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FOREWORD 


This NASA Special Publication is a documentation and discussion 
of a digital computer program used to conduct ray-tracing of 
electromagnetic waves in the magnetosphere. The publication 
consists of three sections. Section I contains a general 
description of the program and its capabilities. Section II 
assumes that the reader is aware of the scope of the program 
and provides all the information necessary for a non-programmer 
to run the program. Section III is a programming manual 
containing extensive information on the program structure. The 
program is designed to operate in the IBM 7094 IBSYS environment. 
It is written in Fortran IV but utilizes a MAP assembler 
subroutine to integrate the differential equations. 

Ray-Tracing techniques are used extensively in the study and 
understanding of propagation of electromagnetic waves in any 
media (neutral/ magneto-ionic , etc.). Combined with experimental 
data, the ray-tracing technique has served as a powerful tool 
in communications research. Data from NASA satellites (ISIS, 
Explorer, RAE, and the like) are better evaluated with the help 
of digital ray-tracing techniques. The present documented 
program has been used by the authors in studying the data from 
the ISIS topside sounder experiments. However, the program is 
applicable to any medium with suitable choice of models. 

Jayaram Ramasastry 

Edward J. Walsh 
April 1969 
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INTRODUCTION 


This report consists of three sections. Each section is 
intended to be an independent treatment of a specified aspect of 
reader interest. 

Section I contains a general description of the program so 
that the reader can gain a general idea of the program's capa- 
bility. 

Section II assumes that the reader is aware of the scope of 
the program and that he intends to use this program in ray trac- 
ing analyses. Section II contains all the detailed information 
necessary for a non-programmer to run the program. 

Section III is a programming manual containing detailed in- 
formation on the program structure. 

This program is designed to operate in the IBM 7094 IBSYS 
environment. The basic program is written in FORTRAN IV but 
utilizes a MAP assembler subroutine to integrate the differential 
equations. 

The program can be considered to consist of five major 

parts : 

(1) An executive routine which governs program flow 

(2) An input section to assess the initial data necessary 
to operate the program. 

(3) An output which prints a history of the path of the 
ray and under option governs the generation of a plot- 
ting tape 

(4) Three mathematical models which characterize the 
electromagnetic properties of the magnetosphere 

(5) An integration routine which evaluates the differen- 
tial equations 

The ray-tracing technique has been used by many people in 
ionospheric and magnetospheric propagation research. Actual ray- 
paths and signal characteristics (like attenuation, path-loss, 
doppler shift, and refraction) in model atmospheres and iono- 
spheres are computed using the ray- tracing program. Because of 
the high accuracy obtainable, the use of a high-speed digital 
computer is preferred to an analog machine for the integration 
of the ray equations. 


The ray-tracing program used in our study is based on the 
Hamilton system of equations as derived in spherical polar co- 
ordinates by Haselgrove and extended by Grossi and Langworthy 
for the investigation of HF and VHF ionospheric propagation. The 
^ r u^ CU ^ ar P ro k^ em concerning our study is the guided propagation 
of high frequency radiowaves along the magnetic field lines of 
the Earth. Some of the echo traces appearing at virtual ranges 
greater than those of the normal vertical incidence echo traces 
on the topside-sounder ionograms have been explained in terms of 
guided propagation of radiowaves along the field-lines. The 
guided propagation along the field-line is made possible by 
field-aligned ionization irregularities (e.g., ducts and shells) 
of suitable . scale sizes and enhancements or depletions. The 
irregularities are assumed to have thickness greater than the 
radio wavelength and hence act as waveguides or "ducts" to trap 
HF energy and produce the long-range echo traces. The guiding of 
rays along field-lines requires irregularities with a certain 
minimum transverse ionization gradient. Propagation between 
magnetic conjugate points or "conjugate ducting" occurs when the 
transverse ionization gradient exceeds the minimum required for 
guidance at the apex of the magnetic line of force. Conjugate 
echoes are recorded by the topside sounder receivers when the 
signal traverses to and from the conjugate reflection points along 
•the magnetic field-line passing through the satellite. Know- 
ing the satellite orbital parameters and the transmitted signal 
parameters, one could conduct ray-tracing utilizing realistic 
models for the magnetosphere. The ray-tracing program yields 
such useful information as the criteria for guidance, the group- 
delays of the trapped signals of various frequencies, and total 
path length traversed. Ray-tracing method is a powerful tool 
when used alongside experimental observations since it gives a 
better insight into the observed results. 
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I. GENERAL PROGRAM DESCRIPTION 


Introduction 

This section provides a general description of the capabili- 
ties of the program. Major emphasis is placed on the program 
input and output and the basic formulation of the differential 
equations . 

The integration technique is discussed only briefly but 
reference is made to Appendix A where the integration package is 
described in detail. 

The three models characterizing the electromagnetic proper- 
ties of the magnetosphere (electron density, collision frequency 
and magnetic field) are discussed briefly in this section and 
their detailed descriptions are presented in Section III. 


Main Program 


The following equations form the basis of the ray tracing . 
program. The first six equations are the Hamilton electromagnetic 
wave equations in spherical coordinates. The remaining five equa- 
tions are additional functions of phase path length. All deriv- 
atives indicated with a dot are with respect to phase path length 
in km. 

The ray is described in terms of position in spherical polar 
coordinates with origin at the center of the Earth and in terms 
of the components of the wave normal y^ , ^2 • Y3 i- n r > 9 / <P 

directions, respectively. Figure 1 describes the ray position 
geometry in spherical coordinates. 



Figure 1. Ray Position in Spherical Coordinates. 
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where : 


y = phase refractive index 
r = geocentric radius of the wave front 
0 = colatitude of the wave front 
<p = longitude of the wave front 

(y r y 2 'Y 3 ) = components of the wave normal in the 
6 r 0 , and t p directions 

G = group path length 

S E ray path length 

D e powerloss due to absorption 

E e coefficient used to compute the doppler 
shift, see Section III 

P p e coefficient used to compute the 
doppler shift, see Section III 

t e eikonal function that has the property 
that the surface t = constant, is the 
geometrical wave front 

a e angle between the wave normal and the 
wave front 

The powerloss calculation, Eq. (11) , is optional and may 
be suppressed if this data is of no interest to the user. 

A complete description of all calculations can be found in 
section III. 


Input Section 

The input routine is written with the aim of providing maxi- 
mum control over the operation of the program and at the same 
time minimize the size of the input deck. The NAMELIST feature 
of FORTRAN IV is used because it gives the user maximum flexi- 
bility and also tends to minimize the careless errors usually 
associated with fixed format input statements. 

The input data is divided into three major parts: 

(1) The data required to describe the initial position, 
direction and characteristics of the ray 
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(2) Optional integration parameters 

(3) Optional limits for the axes of the calcomp 
plots 

The required program inputs are: 


(1) 

Starting position and direction of the electro- 
magnetic wave, which is referred to as a ,: ray n 

(2) 

Integration stops or triggers which 
temporarily or permanently halt the 
procedure 

either 

integration 

(3) 

Propagation frequency, and ray mode 
or extraordinary) 

(ordinary 

(4) 

Output controls, print and plotting 

intervals 

(5) 

Option indicators 


(6) 

Optional input indicators 



All optional parameters are preset in the program. If, 
however, any or all of these parameters must be changed, the 
program allows for easy modification of these nominal values. 

For a complete description of all required and optional 
inputs, see section II. 


Output Section 

All output, with the exception of a list of the case inputs, 
is performed by subroutine OUTPUT. This subroutine not only 
prints a history of the path of the ray, but under option control 
governs the plotting of the data. 

Printed output .- Results are printed whenever one of the 
following situations occur: 

(1) A print time is reached. Print time is a function 
of phase path length. A print time occurs at 

n • Ahp where n = 0 , 1 , 2 , 3 . . . and Ahp is some incre- 
ment in phase path length, measured in kilometers. 

(2) The rate of change of the geocentric radius, r, is 
approximately equal to zero. 
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(3) A stopping condition has been reached. 

(4) A reflection has occurred. 

The printed output consists of a sequence of 3 x 7 matrices 
grouped thirteen per page. The first matrix in the sequence, 
printed at the top of the page, is a heading matrix which identi- 
fies each data element of the data arrays. The heading defini- 
tions are listed below according to row and column. 


Row 

Column 

Title 

Description 

1 

1 

PHASE PATH 

Phase path length, hp, from (r , 0 , <f> ) 
in kilometers. o o o 

1 

2 

RADIUS 

Geocentric radius, r, of the ray position 
in kilometers. 

1 

3 

CO LATITUDE 

Colatitude, 0, of the ray position, in 
degrees. 

1 

4 

LONGITUDE 

Longitude, cp, of the ray position, in 
degrees. 

1 

5 

ABSORPTION 

Absorption loss, D. 

1 

6 

DOPPLER, SP 

Doppler shift 

1 

7 

POWER LOSS 

Pro 

Power loss, , exclusive of absorption 

P T 

2 

1 

GROUP PATH 

Group path length, G, in kilometers 

2 

2 

Y ! 

Vertical component, Y , of the wave 




normal 

2 

3 

Y 2 

Southerly component, Y , of the wave 
normal 

2 

4 

Y 3 

Easterly component, Y , of the wave 
normal 

2 

5 

MU** 2 

2 

Square of the Index of refraction (ju ) 

2 

6 

Y**2 

2 2 2 
Y l +Y 2 + Y 3 

2 

7 

EPSTEIN CD 

Epstein condition 
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Row 

Column 

Title 

Description 

3 

1 

RAY PATH 

Ray path length, s, in kilometers 

3 

2 

POLARIZATION 
-MOD AND ARG 

Modulus and argument of the wave 
polarization term, R. 

3 

4 

DEL MU 

Validity criterion 

3 

5 

N 

Electron density, N, in electrons/cc 

3 

6 

NU 

Collision frequency, v, in collisions/sec. 

3 

7 

GROUP DELAY 
(Gd) 

Group delay in milliseconds. (Group path 
length divided by the velocity of light in 
free space. ) 


An example 

of the printed 

output can be found in Table I. 


Whenever results are printed because of any one of the stopping 
conditions, an appropriate message is printed below the data to 
which the message refers. 

Plotted output .- Each output plot consists of two plots. The 
rectangular plot shows the distance of the ray normal to the field 
line from the near end. The polar plot shows the actual ray path 
with reference to the surface of the Earth. The polar plot gives 
a clear indication of the positions of the conjugate reflection 
points as well as the L value of the field line guiding the rays. 

A complete description of the plotting routine can be found in 
section III. 

Some examples of the plotted output are described in the 
following paragraphs. 


Examples 

Figures 2 through 4 show a few results of the ray tracing. 
They are chosen so as to demonstrate the capabilities of the pro- 
gram. Figure 2 is a single plot which simulates a conjugate echo 
path observed in an Alouette 2 ionogram. Figure 3 is an overlayed 
plot of two rays launched from the same point but with different 
launch angles, DELAO. Figure 4 is an overlayed plot of four rays 
launched from the same point with identical launch angles but at 
different frequencies. 

Parameters shown in the figures are defined as follows: 

RO = the geocentric radius in km of the initial 
signal position 
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NDQ1 = -Iu*b601*t_-i0 


PKFRAC 


= the peak fractional ionization enhancement 
LAMBDA = the colatitude of the field line 

AO = the initial angle the ray makes with the 
local vertical 

BO = the initial angle the ray makes with the 
south vector 

DELAO = the initial angle the ray makes with a line 
drawn parallel to the tangent of the field 
line at RO 

PHI 

& THETA = the magnetic longitude and colatitude of the 
initial signal position 

HO = the scale size of the duct at the base of 
the field line 

FREQ = the frequency of the initial signal 

The ray mode (ordinary or extraordinary)' is also indicated. 

Each figure consists of two plots. The rectangular plot 
shows the distance of the ray normal to the field line versus 
the distance along the field line from the near end. The polar 
plot shows the actual ray path with reference to the surface 
of the Earth. The polar plot gives a clear indication of the 
positions of the conjugate reflection points as well as the L 
value of the field line guiding the rays. 

All cases are for the L = 1.53 field line using a peak elec- 
tron density enhancement of 5.0 percent at the point the ray is 
launched. The rays are all of the extraordinary mode. 

In Figure 2 the ray has an initial geocentric radius of 
9295.40 km. This corresponds to a distance of about 2.75 km 
from the field line at a colatitude of 77.445 degrees. The 
ray starts at about 4,500 km along the field line from its 
base in the northern hemisphere. It passes from north to south 
and is reflected in the southern hemisphere. After reflection 
it travels back into the northern hemisphere where it is again 
reflected. The program was deliberately stopped after, the 
second reflection. 

The sliaht bowing observed in the path of the ray near the 
equator is caused by the larger scale size near the equator. The 
larger scale size causes the enhancement structure to be wider. 
Thus, the ray rides further out from the field line. 
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Figure 3. - Overlayed Plot of Tjwo Rays Launched 
from the Same Point but with Different Launch Angles, DELAO 
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Figure 4B. - Overlayed Plot of Four Rays Launched 




The oscillatory nature of the ray path seen in the rectan- 
gular plot is caused by the ray being launched at a point where 
the electron density gradient is greater than necessary for guid- 
ance. The ray tends to oscillate about a point where the elec- 
tron density gradient is just sufficient to guide the ray paral- 
lel to the field line. A careful launching of the signal can 
minimize these initial oscillations. In this figure the oscil- 
lations were somewhat damped upon reflection. However, this is 
not always the case. For example, in Figure 4 the amplitude of 
the oscillations after reflection was actually increased. 

Figure 3 shows three rays launched from the same point but 
with different launch angles, DELAO. The launch angles corres- 
ponding to plots (a) , (b) and (c) are 0 degrees, - 0.5 degrees 

and + 0.5 degrees, respectively. The oscillatory nature of the 
ray paths in plots (b) and (c) once again demonstrates that the 
rays were launched off the equilibrium position. The rays start- 
ed in the northern hemisphere about 2500 km along the field line 
and propagated along the field line to the southern hemisphere 
down to the reflection level. After reflection, they retraced 
their path and returned to the reflection level in the northern 
hemisphere and were again reflected. The program was stopped at 
this point because two reflections were counted. 

In Figure 4, one can observe the behavior of waves at dif- 
ferent frequencies. Three rays are launched from the same point 
with identical launch angles. The enhancement model is the same 
but the three rays differed in their frequencies. The three fre- 
quencies are 0.95, 1.15, and 1.25 MHz. The ray with a frequency 
of 0.95 MHz has the greatest amplitude of oscillation and the 
ray with the frequency of 1.25 MHz has the minimum initial oscil- 
lations. The rays were launched in the southern hemisphere at a 
colatitude of 77.445 degrees. The rays passed from south to 
north. Only one ray was reflected back into the duct at the con- 
jugate reflection point. The other two escaped upon reflection. 

The inverse proportionality of the oscillation amplitude to fre- 
quency indicates that higher electron density gradients are re- 
quired to guide the rays of higher frequencies. 

Figures 5A, B; 5C , D, and 5E, F correspond to cases where a 
frequency of 1.2 MHz is trapped in enhancement ducts of different 
peak fractional enhancements. The ray launching position and the 
launch angle is the same for all the three cases. However, the 
peak fractional enhancements are 5, 6 and 7 percent for Figures 
5A, 5C, and 5E, respectively. It may be noticed that in Figure 
5A, the peak fractional enhancement of 5 percent is not sufficient 
to contain the ray on its way back from the reflection point in 
the northern hemisphere. Peak fractional enhancements of 6 and 
7 percent as shown in Figures 5C and 5E seem sufficient to trap 
the rays and when the program was terminated after two reflec- 
tions, the rays were still well trapped. 
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Magneto spheric Models 


The following characterizes properties of the magneto- 
sphere : 

(1) The electron density distribution 

(2) The magnetic field distribution 

(3) The electron collision frequency distribution 

The mathematical models for these quantities have been imple- 
mented in the form of subroutines and can be independently modi- 
fied as necessary. 

Magnetic field model.- A dipole model is used for the mag- 
netic field of the Earth. The magnetic field equation which 
defines gyrofrequency , f„, is given by: 

ri 

V r '« = c n (;) 3 [i + 3 cos 2 «] 1/2 

where a = 6378.0, the radius of the earth in km, and r and 6 are 
the geocentric distance and colatitude of any point on the field 
line . 



e 

27rm 


B * 1. OE-6 = 0. 9 
o 


where B 0 is the magnetic field on the surface of the Earth at the 
equator, and e and m are the charge and mass of an electron. A 
value of 0.3142 Gauss is used for B . 


The angle between the magnetic field direction and the 
wave normal is given by 


cosip = 


2Y cos0 + Y sine 
1 2 


E Y i j i 1 + 3cos2 °) 


sin ip = 


2Y cose = Y sine 
2 1 

|2Y cose - Y sine | 
2 1 


1/2 


^l-COS^ 
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TMtrn r 109.050 
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Figure 5A. - Ray Path of 1.2 MHz Signal Trapped in an Enhancement Duct. Peak Fractional 
Enhancement in the Duct at the Initial Ray Position is Equal to 5 Percent 


NORTH 
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Figure 5C. - Ray Path of 1.2 MHz Signal Trapped in an Enhancement . Duct . Peak Fractional 
Enhancement in the Duct at the Initial Ray Position is Equal to 5 Percent 


NORTH 



Figure 5D. - Ray Path of 1.2 MHz Signal Trapped in an Enhancement Duct 
Peak Fractional Enhancement in the Duct at the Initial Ra; 
Position is Equal to 6 Percent 



9295.40 AO - 65.992 
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Figure 5F. - Ray Path of 1.2 MHz Signal Trapped in an Enhancement Duct. 

Peak Fractional Enhancement in the Duct at the Initial Ray 
Position is Equal to 7 Percent 


But, the dipole model is highly idealized and lacks such features 
as local variations in the magnetic field strength. Its advan- 
tage lies in the fact that it saves a lot of computer time and 
for the ray-tracing problem under study, the microscopic features 
of the magnetic field are not necessary.* 

Collision frequency model .- The collision frequency model 
has the following functional form: 


v 


10 


i 


w here v' is computed as shown below. The collision frequency 
profile as a function of altitude consists of three parts that 
are smoothly joined with the aid of a curve fitting program. 

for 6378 < r =$ 6478 km: 


v = 12. 03527 - 0. 07392 x 


where 


x = (r-6378)km. 


for 6478 < r 6853 km: 


= j^a. + C(0, 4>)b. f.(x) 

i= 1 


where 


f x (x) = 1 

f 2 (x) = x 


x = (r-6478)km 


f 3 (x) = X 
f 4 (x) = X 3 

f (x) = cos(0 0 0157x) 
5 


*Those who are interested in using a better magnetic field model 
Should refer to Langworthy (Smith) (1966) for a description of 
the Gaussian spherical harmonic magnetic field model. 
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f (x) = sin(0. 0157X) 
b 



= 5. 0562 

; b = 0. 032512 


a = -3.7482xl0” 2 

Li 

; b 2 = -0.8847xl0" 2 


a = 1.3 864x1 o" 4 

O 

; b = 0.8541xl0" 4 

d 


a . = 1. 4777x1 o" 7 
4 

; b . = -1.5422xl0“ 7 
4 


a r = -0.48192 
5 

; b = 0.01470 
5 


= -0.27021 

6 

; b = 0.65037 
6 


c(0> <f>) = c x e 2 + c 2 0 + c 3 

+ (d^# 2 + d 2 0 + d g ) coS(f> 

where 

C = -0.35818 

; & 1 = -0.17828 


C 2 = 1.1250 

; d 2 = 0.55997 


C 3 = -0.88344 

; d g = 0.56028 


0 is the colatitude and eft is the longitude in degrees. Symbol 
<J> = 0 corresponds to local noon. 

For r > 6853 km, 

v = 2.3653-0. 0030266X + (0.3195-0. 0000536x) C(6,<p) 


where 

x = (r-6853) 

Electron density model .- The model for the electron density 
distribution incorporates those features that have a bearing on 
the high-frequency ducting problem. Thus, an ionization ir- 
regularity (duct) model is superimposed on the normal radial 
distribution of electron density. The duct structure is aligned 
along the magnetic field. 

A complete description of the electron density model can. be 
found in section III. 
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Integration Package 


Because of the high accuracy obtainable, the use of a digi- 
tal computer is preferred to an analog machine for the inte- 
gration of the ray equations. The method used solves the first 
n of a set, N, of first-order differential equations simultane- 
ously. The Adams-Moulton open and/or open and closed formulas 
are used. A Runge-Kutta fourth-order integrator is used as a 
starting routine to generate the necessary differences. Pro- 
vision is made for interrupting the integration process at speci- 
fic values of either the independent or dependent variables. The 
order of differences used in the Adams-Moulton mode is less than 
or equal to nine. 

A complete description of the integration package is found 
in Appendix A. 
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II. USER'S MANUAL 


Introduction 


tion Ind S rrnf!Ti® Manual is intei >ded to supply all the informa- 
tion and guidelines necessary for a non-programmer to operate 
the ray-tracing program. p ate 


Inputs 

In order to facilitate using this program, 
feature of FORTRAN IV is used. This input mode 
advantages over a fixed format statement: 


the NAMELIST 
has two distinct 


(1) It minimizes keypunching errors 

(2) When stacking cases, it minimizes the size 
of the input deck 


NAMELIST allows for a free format input stream; the data 
can be placed anywhere on the card and spaced in any convenient 
The basic *- ul f s the us er must keep in mind when con- 
tructmg an input deck are outlined in section III. For a de- 
tailed description of the NAMELIST feature see IBM form C28-6390. 


IBM 7090-709 4 IBSYS operating system, version 13, FORTRAN IV 
language . - It has been found that when stacking data sets for run- 
ning as many cases as possible, there is usually some similarity 
among sucessive sets. NAMELIST allows the user to take maximum 
advantage of any similarity. 

The first data set in the stack must contain all data re- 
quired by the program for normal operation. Each successive 
data set, however, must specify only those parameters whose values 
differ from immediately preceding case. 


We consider the input data to be of three distinct types: 

(1) Data which is required for the normal operation 
of the program and must be supplied by the user 

(2) Optional integration control parameters 

(3) Optional plotting limits 

Detailed descriptions of these three data types can be found in 
the following paragraphs. 

Input dictionaries .- The input is under control of three 
dictionaries. The first dictionary, XNAME1, refers to the re- 
quired program inputs. The remaining dictionaries, XNAME3 and 
XNAME5 , refer to the optional integration and plot parameters, 
respectively. 
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Required input - XNAME 1 . - Figure 6 shows a sample input deck 
for the required inputs. The dictionary for these inputs consists 
of five names: D, LAST, NTITLE , NCONST, and LIMITS. The first 
name, D, is a 28-word array containing all the data about the 
initial position of the ray, direction and ray characteristics. 

The four remaining names are single word indicators. An ex- 
planation of each of the required inputs follows. 

D array.- Table II contains an item by item description of 
each of the parameters in the D array. The subscript shown in 
the left-hand column indicates the position in the array. The 
subscript of the first item must be punched and then the remain- 
ing items will be stored sequentially. Figure 7 is a graphical 
illustration of the initial ray-position input parameters speci- 
fied by D (1) through D(6). 


TABLE II 
D ARRAY 


D 

DESCRIPTION 

1 

Initial geocentric radius, r , in kilometers 

2 

Initial colatitude, 0 , in degrees. (0^64 180) 

3 

Initial longitude, <|> , in degrees. (0 >< M 360) 

4 

Initial ray angle, A , with respect to the local verti- 
cal, in degrees. (0° £ A q ^ 180°) 

NOTE 

A q is no longer a program input. 

A q is calculated from AA q . How- 
ever, the location D(4) is still 
reserved for A q . 

5 

Initial ray angle, B , with respect to the south vector, 
in degrees. (0° < B q 360°) 

6 

Initial ray angle, AA , with respect to the tangent of 
the field line in degrees.' (0 « AA q 180 ) 

7 

Maximum allowable geocentric radius, R max > in kilometers 
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TABLE II.- Continued 



Description 


Minimum allowable geocentric radius, R min > in kilometers 

Maximum allowable colatitude, 0 m3 , in degrees. 

(0° * 0 max 4 180 } 

Minimum allowable coatitude, 6 . , in deqrees. 

(0° 4 9 min 4 18 ° 0) min 

Maximum allowable longitude, Vax ' ln degrees, 
max " ' 

Minimum allowable longitude, d) . , in degrees. 

(0° « <f) min 360°) min 

Print interval in kilometers. Printed output is keyed 
to phase path length 

Plot interval in kilometers. Plotted output is keyed to 
phase path length. It has nothing to do with the Print 
Interval. It is only a control of the spacing of the 
points to be plotted. 

Nominal integration' step size, in kilometers. Integra- 
tions are with respect to phase path length 

Wave frequency, F, in megacycles per second 

Ray type indicator: 

Set equal to 1 for ordinary type, or 
Set equal to -1 for extraordinary type 

Reflection indicator: 

Set equal to 2 if reflection is desired, 0 otherwise. 
Ordinarily, reflection occurs whenever the rate of 
change of geocentric radius goes to zero 

Power loss calculation indicator: 

Set equal to 1 if computation is desired, 0 otherwise 
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TABLE II.- Concluded 


D ARRAY 


D Description 


20 Plot option indicator: 

Set equal to 1 if plotting is desired, 0 otherwise 

21 Plot overlay option indicator: 

Set equal to 1 if overlaying of successive plots is 
desired, 0 otherwise. Ignored if NPLOT = 0 

22 Automatic positioning of ray at point of maximum elec- 
tron density gradient indicator: 

Set equal to 1 if automatic positioning is 
desired, 0 otherwise 

23 Scale size of the ionization irregularity at the base 
of the field line H 

o 

24 Peak fractional enhancement at (r Q , 6 0 , 4> 0 ) where r Q 
is the initial geocentric radius, 0 O is the initial 
colatitude and <j) 0 is the initial longitude 

25 Colatitude, A, in degrees of the field-line passing 
through the initial ray position. 

26,27,28 Not used 
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Figure 6. - Sample Deck Setup 


EARTH'S 

SURFACE 


RAY-POSITION 

( r o» 


LOCAL 

VERTICAL 


AA 0 

RAY- 

Zs VECTOR 


MAGNETIC 

FIELD 


y 2 

(SOUTH) 


0 ^ 

EARTH'S 

CENTER 

Figure 7. - Initial Ray - Position Inputs 




of XNAMEl are four indicators described in Table III. 


TABLE III 

REQUIRED INPUT INDICATORS 


Name 

Description 

LAST 

Last Case Indicator: 


The program assumes LAST =0, so the 
last case to be executed must explic- 
itly set this parameter to 1 

NTITLE 

Optional Title Information: 


NTITLE preset to 0 . If some identifi- 
cation is to be printed, set NTITLE = 1 

NCONST 

Optional Input Indicator, Integration Para- 
meters : 


The program assumes NCONST =0. If 
the integration parameters are to be 
changed, set NCONST = 1 

LIMITS 

Optional Input Indicator, Plotting Limits: 


The program assumes LIMITS =0. If 
the plotting limits are to be changed 
set LIMITS = 1 


Optional input .- The parameters described here are preset by 
the program. However, the user can replace any subset of these 
parameters by setting the values of NCONST and/or LIMITS equal 
to 1 and then punching a set of data cards using the dictionary 
XNAME3 and/or XNAME5. These cards must follow the statement 
card if NTITLE = 1 or immediately following the required in- 
puts if NTITLE = 0. 

Integration parameters - XNAME 3 . - The optional integration 
parameters are under control of the dictionary XNAME3. The 
dictionary contains seven names representing seven data items. 
Table IV defines each parameter and gives its nominal value. 
Appendix A contains a description of the integration routine. 
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TABLE IV (XNAME3) 


Integration Parameters 


Name 

Nominal 

Value 

Description 

ORDER 

1.0 

Precision option. When set to 0 the 
integration is carried out in partial 
double precision, set to 1.0 integra- 
tion is performed in full double pre- 
cision 

EUBAR 

1.0E-5 

Maximum integration error. When the in- 
tegration error is greater or equal to 
this value the integration interval is 
halved 

ELBAR 

1.0E-7 

Minimum integration error. When the in- 
tegration error is less than or equal to 
this parameter the integration step is 
doubled 

YCLOW 

1.0E.7 

Minimum value of the dependent variable. 
If the dependent variable is less than 
or equal to this parameter no halving 
will take place 

HMAXT 

1.0E+4 

Maximum integration step. The integra- 
tion step is not permitted to exceed this 
value 

HMINT 

1.0E-7 

Minimum integration step. The integra- 
tion step is not allowed to be smaller 
than this input 

KD 

4 

Type of integration (must always equal 4) 


Plot limits (XNAME5 ) . - The optional plot limits are under 
the control of dictionary XNAME5. This dictionary contains 
eight names. Table V defines each limit and its nominal value. 
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TABLE V (XNAME5) 


PLOT LIMITS 


Name 

Nominal 

Value 

Description 

XMINO 

O 

o 

Minimum value of the X-axis 
for the rectangular plot. 

XMAXO 

1.20E4 

Maximum value of the X-axis 
for the rectangular plot. 

YMINO 

o 

00 
I — 1 
1 

Minimum value of the Y-axis 
for the rectangular plot. 

YMAXO 

+2.0 

Maximum value of the Y-axis 
for the rectangular plot. 

XMIN1 

O 

O 

Minimum value of the X-axis 
for the Polar Plot. 

XMAXl 

2. 0E4 

Maximum value of the X-axis 
for the Polar Plot. 

YMINl 

o 

o 

Minimum value of the Y-axis 
for the Polar Plot. 

YMAX1 

1. 0E4 

Maximum value of the Y-axis 
for the Polar Plot. 


Input conventions .- The following is a listing of a few sim- 
ple rules the user must keep in mind when constructing an input 
deck : 


1. Column 1 of each input card is not to be used. Any infor- 
mation punched in this column will be ignored. 

2. The first character encountered by the input system must 
be a dollar sign ($) punched in column 2 of the first card 
of each group of inputs. 

3. All data must be separated by commas. The input system 
assumes that any single input quantity must lie between 
an equal sign and a comma, or between two commas. 

4. All numerical data can be entered in either a fixed decimal 
or exponential format. Whenever a decimal point does not 
appear within a number, in the case of the fixed decimal 
format, the system assumes that it is in front of the sep- 
arator and all trailing blanks (i.e., no punches) will be 
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converted to zeroes. In the case of the exponential 
format the decimal point is assumed to lie to the right 
of the least significant digit. As an illustration, the 
following list of numbers will have identical binary 
representations in core storage: 

2500.0 

25bb fixed decimal (b-blank, no punch) 

25E+2 

2.5E+3 exponential format 

Note that in the above illustration the exponent is 
separated from the significant part of the number by 
the letter E. 

5. The NAME LIST system allows for continuation cards. If 
continuation cards are used then the user must be care- 
ful in placing the input cards in proper order. When 
using continuation cards the name of the table, (array) 
need not be repeated. The only limitations are that 
column 1 is not to be used and each card must end with 
a comma . 

This program allows the user to execute several cases 
yithin a single run. This is done by first setting the 
indicator, NA, equal to zero within the first set of 
input and then for each of the following cases simply 
punch the dictionary name with the proper system sym- 
bols and then specify only those inputs which are to 
differ from the immediately preceding case. The inputs 
for each case are then placed behind each other, no sep- 
arator cards are required, and the resulting input deck 
is placed behind the binary deck as described earlier. 

6. The end of an input string is indicated by a dollar 
sign ($). The input system will continue reading until 
the second dollar sign is encountered. 


Operating Instructions 

Deck setup .- The following is a list of the control cards 
necessary to compile and execute the program in source form: 

CC CC CC 

1 8 16 

$ JOB D3121A 1 ERC J. Ramasastry - TRACE 
$EXECUTE IBJOB 

$ IB JOB FIOCS 


35 


$ IBFTC MAIN 


$ IBFTC INPU 


$ IBFTC OUTPU 


$ IBFTC FIEL 


$ IBFTC DENS 


$ IBFTC COL 


$ IBFTC CALCO 


$ IBFTC SCA 


$ IBFTC FORC 


$ IBFTC POLA 


$ IBFTC POWER 


$ IBFTC MARS 


$ IBFTC SMK2 


MAIN PROGRAM 


SUBROUTINE INPUT 


SUBROUTINE OUTPUT 


SUBROUTINE FIELD 


SUBROUTINE DENSE 


SUBROUTINE COLL 


SUBROUTINE CALC04 


SUBROUTINE PRAM 


SUBROUTINE FORCE 


SUBROUTINE POLAR 


SUBROUTINE POWERL 


SUBROUTINE MARK 


SUBROUTINE SMARK 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


(SOURCE) 


$IBFTC CST 


$IBFTC MIN 


$IBFTC TO 


SUBROUTINE CSINT 


SUBROUTINE MINV 


SUBROUTINE TOR 


$DATA 


DATA DECK 
' END OF FILE CARD 


(SOURCE) 

(SOURCE) 

(SOURCE) 


I/O assignments .- The program uses the standard ERC-IBM 
7094-11 assignments. The FORTRAN IV logical units are 5 and 6, 
respectively. If a plot is to be generated a tape must be mount- 
ed on SYSUT7 . " 

Abnormal termination .- The program recognizes three types 
of abnormal terminations : 

(1) An error in the integration routine. The message, 

ERROR RETURN FROM MARK is printed. The program will 
continue on to the next case. 

(2) y 2 ■$ 0 

The message, THE VALUE OF EMUS IS NEGATIVE is printed. 
The program will continue on to the next case.' 

(3) 1 - x ^ 0 

The message, THE VALUE OF TERM IS NEGATIVE is printed. 
The program will continue on to the next case. 
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III. PROGRAMMER'S MANUAL 


Introduction 

This Programmer's Manual contains all the technical infor- 
mation required by a programmer whose task it may be to modify 
the program. 

Each subroutine is described separately. Each description 
generally consists of: 

(1) The mathematical formulation of the computations 

(2) A flow chart of the routine 

(3) A dictionary of the FORTRAN variables calculated in 
the routine 

Following the subroutine descriptions the labeled COMMON 
statements are correlated with the subroutines in which they are 
used. 


Source statement lists of all the subroutines are provided 
in Appendix B. 


Main Program 

Description . - The main program, MAIN, serves as an executive 
routine. (See flow chart in Figures 8A to 8H.) The functions 
of MAIN are : 

(1) Initialization of the integration routine 

(2) Computation of the derivative box 

(3) Handling of trigger stops 

The integration routine has a return indicator, NRTN, which, 
when tested, determines the order of calculations. 

Upon return from the integration routine, the variable 
NRTN is set to an integer value from 1 to 5. NRTN is used to 
control program flow. Table VI indicates the action taken 
when NRTN assumes a particular value. 
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Figure 8B. - Flow Chart of Main Program 
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COMPUTE ELECTRON 
DENSITY (EN) , 
AND DERIVATIVES 


COMPUTE 
INTERMEDIATE 
TERMS USED 
BY SUBROUTINE 
FIELD 


DENSE 



COMPUTE ELECTRON 
COLLISION 
FREQUENCY, (GNU) 


COLL 


Z = (1.59E5) 

r 


' COMPUTE: \ 
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DERIVATIVES 


Figure 8C. - Flow Chart of Main Program 
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Flow Chart of Main Program 
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Figure 8E. 


Flow Chart of Main Program 
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Figure 8G. - Flow Chart of Main Program 
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Figure 8H. - Flow Chart of Main Program 



TABLE VI 


ACTION TAKEN WITH REFERENCE TO RETURN INDICATOR 



Trigger stops are used to interrupt the integration routine 
at predetermined points. This program uses ten triggers. If a 
trigger stop occurs (i.e., NRTN = 4), the trigger indicator, 

NTRG , is set to an integer value from 1 to 10. The value as- 
sumed by NTRG corresponds to the trigger activated by the in- 
tegration routine, MARK. 

Table VII correlates the value of NTRG, the trigger acti- 
vated and the sequence of actions taken. 
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TABLE VII 


SEQUENCE OF ACTIONS TAKEN WITH 
REFERENCE TO TRIGGER INDICATOR 


NTRG 

1 

Print interval: call OUTPUT to print 

data, update trigger, return control 
to MARK. 

2 

The geocentric radius, r, has reached 
R max : call OUTPUT to print data, if 

NPLOT > 0 plot data, print end of case 
message, call INPUT to process next 
case. 

3 

The geocentric radius, r, has reached 
R m in : if reflection indicator, JTEST, 

does not equal two, the procedure is 
the same as for NTRG = 2, if JTEST = 2, 
the ray is reflected from a plane per- 
pendicular to r and control is return- 
ed to the integration routine. 

4 

0 > 0 max : the Procedure is the same 

as for NTRG = 2. 

5 

0 ^ 9min : the procedure is the same 

as for NTRG = 2. 

6 

<(> > <|) max : the procedure is the same 

as for NTRG = 2. 

7 

<J> < <j>min : the procedure is the same 

as for NTRG = 2 . 

8 

Phase path, h p , has reached a value in- 
dicated by TEST J : the stop on TESTJ 

has two functions. The first is to 
delay the start of the power loss cal- 
culation until the ray has traveled 
enough to be considered in the far 
field. It also turns on the trigger 
for a reflection stop when the ray has 
gone a specified distance past the re- 
flection. 




TABLE VII.- Concluded 


SEQUENCE OF ACTIONS TAKEN WITH 
REFERENCE TO TRIGGER INDICATOR 


NTRG 

9 

r = 0: call OUTPUT to print data, 

and store current values of r, 0, 
and 4> in plotting arrays, print 
value of r, return control to MARK. 

10 

Plot interval: store current values 

of r, 0, and <f> in plotting arrays, 
update trigger, return control to 
MARK. 


Differential equations evaluated in the MAIN PROGRAM . - The 
following equations form a closed system for tracing the rays. 

The ray is described in terms of position in spherical coordinates 
with the origin at the center of the Earth and in terms of the 
components of the wave normal, in the r, 0, 4> di- 

rections, respectively. 


#1 


(' 

jU \ 


r = — ( Y -a 

° 1 IN M 9Y 


; r(o) = r 


#2 e = — n 


( Y - u 
2 l 2N M 9Y 
TfX \ 2 


; 0(o) = 0 


#3 0 = 


r— ( 

sm0 \ 


2 • „ 
ryj sm0 


V 9k 

Y ~ix — 
3N M 9Y„ 


; 0(o) = 0, 


« Y 1 = t ft + » Y 2 + *Y 3 »toe 


• Y i<°> = Y io 
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#5 



& 

90 


“ i*Y + $Y rcos0 

^ O 


Y 2<°> * Y 20 


#6 


Y = — f — 
3 rsin0 


10 u . 

; ^ - rY 3 sln « - grease 


Y 3<°) ■ Y 30 


where Y in are components of the normalized Y-vector, i.e.. 




i= 1 

Y is normalized by subroutine FORCE. The symbols r, 0, <p , Y lf 
y 2 ' y 3 are derivatives with respect to phase path length, 
hp (the independent variable), in kilometers. The symbol y is 
the phase refractive index. 

The simplified Appleton-Hartree expression for the phase 
refractive index is used. It is given by 


2 




= 1 - 


X 


1 - M + 



and the collisions are neglected. 

The appropriate algebraic sign (+) in the denominator is 
chosen according to the indicator supplied by the user via the 
required input D(17). The plus (+) sign refers to ordinary rays 
and the minus (-) sign refers to extra-ordinary rays. 

2 

v _ e N _ N N 

2 o p 

€ o m(27Tf) 12400. f 1 

N = electron density in electrons per cc 

f = frequency of propagation in megacycles 
per second 

e = charge of an electron 
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m = mass of an electron 


= permittivity of free space 
e m(27rf) 2 

C x = ~ 2 — ^ = 12400. i 


M = 



1 - X 


T 


f H sin \p 

f 


f nos ih 



f = gyrofrequency in megacycles per sec 
H 


ip = angle between the magnetic field and the 
wave normal 


The problem is now reduced to one of finding N, fR and ip 
and their derivatives with respect to r, 0, 4> , Y^, Y 2 and Y 3 . 

The initial coordinates of the ray-position (r 0 , 9o' 4*o) and 
the initial ray direction (Y i0 , Y 20 , Y 30 ) and the wave- 
frequency, f, should be known. Thus, for a given electron densi- 
ty (N) and magnetic field (fR and \p) , one can trace rays of any 
frequency from any point. This is subject to the constraint 
that > 0 . The local values of electron density, N(r, 0, 4>) , 

and the magnetic field, fR(r, ® > $) and ^ ( r ' $ > ^1 ' Y 2 • ^3^ 

are derived from the respective models. 


For Hamilton's equations, one must evaluate 9y/9r, 9y/90, 
9y/9<f> , 9y/9Yi, 9y/9Y 2 and 9y/9Y3- 


9y _ 1 

9r 2yM 

K 


_L M + N 

c i 8r c i m r 


9M 


9r 


+ 1 

+ i 

9f 

+ COS 

f 2 

* != i. 

cos ijr 

(] 

R 

\ 3r f H 

9r 

f 2 

9r j 

f 
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9/u _ 1 

90 


R 


1 

9N , 

N f 

9M 1 / 0M 


90 + 

c i m r L 

.‘M 1 5 ( M ^ + 

4 

9cos $ } 

ill 



y! 


9f_ 


.2 90 


ik _ 

1 

( 1 9N N 

r 9M i / 

d(p 

2 wMr - 

1 °1 ^ C ! M r 

[- H 1 5 ( 


Y? 


9f. 


*™*| i ?)]| 


s- k K [-% (“?*-• ?«*)]} 


9^ _ 1 


9y, 


2 2/kM 


_§tL _ 1_ j._N 

9y„ ~ 2 

3 2 M M^ 


_ ( _N f 9M , 1 / 

i I C 1 L 8y 2 - 5 1 M «* 


9M 
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f 2 

8Y i/ 

4 

9C0S# 

f 2 

% 

f 2 

_H 

dcosil) \ 




9M 

9r 


9f 


2f 2 (l-x) 


V in2 ^ + f tsi#^ 
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H . 2, 9N 
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9M 

90 




2 . 9f H 


2f„ sin ip — + ij.simp 


,2 . , 9sini b 


H 


90 H 90 


f 2 

H . 2 9N 
Sm * 90 

1 - X 


9M 


2f (1-X) 


f 2 

H . 2, 9N 

" 5 " sm ^ 

, . 9sin$ , 1 

l2£ H S ‘ n '* 8* + V*"* 8* 


. 2 M H j2 

5in ib -* — + i„ 


-2M _ — L_ Ji simp MM 

9Y (1-X) f 2 V 9Y 1 


_ _1_ Ji Bln * Mat 

9 Y _ 1-X f 2 * 9Y 2 


9M 

9Y„ 


1-X 


f 2 

H . . dsmip 

7 S1# «, 


2e 

9f 


^uM 


Y T Y t ! / Y t Y L 

1-X R l 1-X *T L 


R 


where 




R = if M 2 + Y^ 


The derivatives of 

The derivatives of 
(subroutine) . 


m r = 1 - M + R 

N are calculated in DENSE (subroutine) . 

\p and f are calculated in FIELD 
H 
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Other differential equations evaluated .- Group path length, 
ray path length, power loss due to absorption, doppler shift, 
and power loss along the ray's path are the other differential 
equations evaluated and are discussed as follows: 

(1) Group Path Length 


G 


L 


1 


f 9u 
y 9f 


G l (o) = 0 


G is measured in km 
Li 


Group delay is determined by 


G 


D 



L 


where 


c = 3.0 x 10 5 km/sec = velocity of light 


(2) Ray Path Length 

1 


S = 


y cos a 


S (o) = 0 


S is measured in km 


a is the angle between the wave normal and 
and the ray direction 


(3) Power Loss Due To Absorption 

D = - 2 - D D(o) = 1.0 

y 

„ _ 2irfk 
K - 

where 

D = power loss due to absorption 
k = index of absorption 
K = coefficient of absorption 
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k = 


BX 

3u(A 2 + B 2 ) 


X = 


N e 

c 

e o m(2 Trf)' 


\ "f sin%(l-X) 

A - i - g 2 - 

(1-X) + Z 


U + (U 2 +V 2 ) 


1/2 


. \l/2 
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I Z f sin 2 * 

± 

(1-X) 2 + z 2 


V 


(x 

2 2 1//2 

J 1/2 

)2 

U + (U + V ) 



- 
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12 2 

2 

"l „ 2 . 2 , “ 

2 y sin 


- Z y sin ip 

2 2 
(1-X) + Z 


2 2 
(1-X) + Z 


2 


+ y cos$ 


V = 2 


12 2 12 2 

2 y sin ^(i-x) 2 y sin 


2 2 
(1-X) + Z 


(1-X) 2 + z 2 
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where 



Z = v = collision frequency in 

collisions per sec 


v is given by the collision frequency model. The collision fre- 
quency model is described in subroutine COLL. 

(4) Doppler Shift 

Af = - ^ jucosa ds ; Af(o) = 0 

r 

o 

Af is the Doppler-Shif t in MHz 
c is the velocity of light 
f is the wave frequency om MHz 
y is the phase refractive index 

(5) Power Loss Along the Path of the Ray 
The total power loss in dB is: 


10 log 10 


2 r Mr,- 1 

c o Final 

— jx — x — EXP(E ) x EXP(E ) 

4rrf F T 1 2 



o 
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where 


ds 


dh 


y cos a 

s = ray path length 



= phase path length 


= isotropic radiator aperture loss 


= near field loss due to refraction 
and inverse square law spreading 



far field loss due to refraction 
and inverse square law spreading 


EXP - 


“abs ds = power loss due to absorption 


a ABS related to the absorption coefficient, K, of the 
wave by the expression 


“abs ” K OOS “ • 


Then the terra 



is equivalent to D, the power loss due to absorption. 

Power loss in the near field.- The power losses in the near 
field is calculated in MAIN. It is evaluated only once, when 
the ray has traveled far enouqh to be considered in the far 
field. 
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Pro/ p T gives the power loss in the near field from s 0 to 
s t due to refraction and inverse square law spreading. The equa 
tion given below is for a dipole antenna. 

r and y are shown in the following Figure 9. 



oo 

C . (P) = . f S£^ dx 

'P 

C ± (P) = [ 2lsa 

Jo 

hen 




where 


T = 7TjU 


2 2 2 

(rM = r + r - 2r rcosy 
o o 


cosy = sin(0 o )sin(0)cos(0 o -0) - cos((? o )cos(0) 


y ' = 7T - COS 


-1 


r 2 2 4 2 “1 

( r ') + (r Q ) -r 

2r»r 


A = . 5772 + ln(T) - C.(T) 


\ sin(T) [ S.(2T) - 2S.(T)' 


+ \ cos(T) [C.(2T) - 2C.(T)] 
+ | cos(T) [ . 5772 + ln(| T)1 


Program Notes.- The program variable PROPT is actually 


ro 


9 19 

47rf /a x 10 


where 


c = the velocity of light in km/sec 

5 

= 3 x 10 km/sec 

Power loss in the far field . - The term 

s„ 


^ Final 


v 2 


ds 
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where £ is the eikonal function which has the property that 
the surface £ = constant is the geometrical wave— front and 
gives the power loss in the far field due to refraction and in- 
verse square spreading. These are field losses calculated in 
POWERL routine but integrated under control of the main program. 

Polarization Term .- The expression R = Err/E^r, a complex 
number, is the polarization term. and E^ the compo- 

nents of the electric vector of the wave along axes (x,y) in 
the wave front, y being in and x perpendicular to the plane 
containing the magnetic field of the Earth. 


|R | = modulus of R = 



JY 4 + A + Y 2 A 1/2 
4 T — T 



$ = argument of R = tan 


-if - AU2 [ cos | < 1 ' X >- Z sin |] 


2 ZY T- + A 


1/2 


sin|(l- 


X) + Z cos 


J / *1 Tn 2 2 

d = (1-X) + Z 


A = 


-Y 4 + Y 2 d 
4 T L 


-i 2 


-,2 


1/2 


2 Z Y t (1-X) 

Li 


where 


8 = tan 


-1 


-2(1 -X)Z Y 
4 T L 


Evaluation of the polarization term is described under sub 
routine POLAR. 

Dictionary of Major FORTRAN Names .- Table VIII contains a 
dictionary of major FORTRAN names for the MAIN program. 
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TABLE VIII 


FORTRAN 


NAME 

DEFINITION 

PRNTO 

Value of PRINT trigger 

PLOTT 

Value of PLOT trigger 

MORDGR 

Order of differences carried 
in integration routine 

MUFLAG 

Pass indicator subroutine 
POWERL. Equals zero for 
the first pass and 1 for 
the rest of the run 

DMDR 

3M 

3r 

DMDT 

3M 

3-0 

DMDP 

3M 
3 <f> 

DMDY1 

3M 

3Yi 

DMDY2 

3M 

3Y 2 

DMDY3 

3M 

DMDF 

3M 

3F 

DMUDR 

3y 

3r 

DMUDT 

3y 

30 

DMUDP 

9y 

3<f> 

DMUDY1 

3y 
3 Y x 

DMUDY2 

3y 

3Y 2 

DMUDY3 

sy 

3Y 3 
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TABLE VIII.- Concluded 


FORTRAN 



NAME 

DEFINITION 

DMUDF 


3f 

CAPPA 

K, absorption coefficient 

EPSTIN 

Epstein condition 



3M 

DMDSI 



DMUDSI 


dip 

YD 

Array of values of differ- 


ential 

equations 

DEN0M 

A 1 1 


C0SPSI 

c0s y | 


RP 

r' \ 

see near field 
power losses 

PSI 

Y' j 


PR0PT 

P ro ' 

P 



T 

' 


Subroutine POWERL 


Description . - POWERL computes the power loss in the far 
field due to refraction and inverse square law spreading. 


POWERL evaluates 


P F ■ 


v 2 ; 


2 

fx cosa 

where £ is the eikonal function which has the property that 
the surface £ = constant is the geometrical wave front. 
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4 > 


cos 6 80 

sin0 3r ' 


The derivatives of r' with respect to r, 0, <J> are: 
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To find the derivatives of r, 0, <f> with respect to r', we 
differentiate the equations defining r', 0', 4> * with respect to 
r'. In general, we have 
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- r sin0 sin<£> 7^7 

We then solve these three simultaneous equations for 
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Now we must find the derivatives of y with respect to the 
r, 0 , (J> system. The first derivatives of y are given in MAIN 
PROGRAM. It remains to find the second derivatives. 
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3 2 M 

dtydcf) 




£ 


sin0 


3 2 M 


Sy^r 


_1_ 

2M 


_ 2 MM. 

dr dy 


3 2 U 

8y 2 36 


2m 


a 2 c 


_ 2 & JtL 

30 3y_ 


9 2 u 
3y«,3 0 


2M 


a 2 c 


3f 

f 2 

b _H 

^ 0 * 1 j 

& 30 

2 sin^ 
f 

3cos 0 \ 

N 

90 / 

2mC 1 A 

/ 3A 

3N 2N 

K 

3r A 


30 


3A 3A 


3r 


+ 


3A 

0yo 


/ 3N _ 2N 3A 
\ 30 A 30 


)* 


1 3A / 3N _ 2N 3A \ 

2 n 3y q V 30 A 30 / 


N 

3 2 A 

A 


N 

9 2 A 

A 2c i 

3y 2 30 

N 

3 2 A 

. 2 ^ 
AC i 

3y 3 90 


_ 2 &L JE£_ 
30 3y 
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where 


and 


A = 


i - m ± 


f 2 

M 2 + cos% 

f 


B = 


V 


f 2 

M 2 + -? eos% 
f 


M= 2 T 


f 2 2 
II sin ip 

C 


c=1 _ JI 

C 1 


9 2 M 

dipdr 


2f 




8f H ^ ^ . 9sin?/> 

4 f R sin ip cosip — + 2 fjj cosj/i — "-*■ 


2 . , aces# 2 f l sin ^ 003 ^ If 

+ 2 4 sm * te + 


( c i [ 1 ' ^7] ) 


9 2 M 

dlpdQ 


^ R 


8f H 2 

4f sirup co sip + 2f__ 

H V * 90 


. 9sin^ 

h cos ^ ~^r 


+ 2 


4 sin ^ 


2 9N 

. , soogi , 2£ h sini ^ 00B » w 


90 
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9 2 M 



2f H ^ 9y 


9sin $ 

9r 



8sin# 

9r 


asm# 

*1 


+ 



a 2 sin# 
ar3y x 


+ 



sin# 


8sin# 

*1 


8N / 
3r I 


2 

a m _ i 
"<*2 = 



sin if> 


8sin# 

*2 



3sin# 

% 


9 sin# 
90 



a 2 sin# , 

aaay 2 



sin# 


8sin# 

^2 




£a_ 

drdy 


-OL + -±- ( M ~ + Y cos ip 

drdy b 3 \ 8r f 2 


9 f H 
cos ip — 


+ f . 


dcqstp 


H Or 



t 


9M ~H , 9cos ; p 

M _ + - COS* ^ 


o 2 f 9f 

1 ) m 9M + J_M + _JH COSi/) 

— B J 9r 9y 1 9r9y 1 f 2 d Y 1 ^ 


*H I 9cos j 9cos^ + . 

9r 9y x 9^ 


9 2 A 


9 2 M 


= " + [ M 1ST + "f" C ° S ^ 

909y 2 999y 2 B 3 \ d 9 ? 


COSIp 


90 


+ f 


9cosift 
H 99 



M 


f 2 

9M _H 

^2 f 2 


COS ip 


9cosi/) \ 

** j 


9 2f 9f 

+ 1 ) 8M M + M -£*- + -# COS* ^ -J5 

i B ) 89 0y, 8S3y 2 f ^ ® 


a 


H 


( 222*4 22224 + COSi) , ^2t 

V 90 9y 909y 2 
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9 2 A _ _ 9 2 M - 1 ( tlt 8M f H 

d<pdy 3 ~ d<pdy 3 + fi 3 ^ M d<p f 2 




1 

B 


9M 

d<t> 


9M u. ™ 9 2 M 

+ M — 

9y 3 909y 3 


+ 


2f H 

— cos ip 


dcosjj) 

% 


+ 




dcosip 

% 


+ cosip 


2 

9 cos ip 

9 09y, 


The derivatives of fjj, cos \p, and sin \p' depend on the 
particular magnetic field model used. 

The power loss as computed above (and ray-tracing in gen- 
eral) is valid only when the following conditions are fulfilled: 


J_ Ym 

kM- At 


<< 1 where k = 


2 £ 

X 


and 


k|i 


VPL 

PL 


« 1 


These conditions are violated in rapidly varying media or 
when a caustic is encountered. Figure 10 is a flow chart of the 
subroutine POWERL. 


Dictionary of major FORTRAN names. -Table IX contains a 
dictionary of major FORTRAN names, subroutine POWERL. 
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Figure 10A. 


Flow Chart of Subroutine POWERL 






0 


Figure 10B. 


Flow Chart of Subroutine POWERL 




Figure IOC 


Flow Chart of Subroutine POWERL 



FORTRAN 

Name 


Definition 


DRPDR 


DRPDT 


3r ' 
3r 


3r 1 
39 

9r' 


DRPDP 


DSITP 


COTP 


COPP 


DRDRP 


DTDRP 


DPDRP 


30 


COS (0») 


COS(9 ') 


sin(9) 

3r 
3r ' 


39 
3r ' 


30 


D2MY1R 


D2MY2T 


D2MY3P 


3r' 

3 2 M 
85 ^ 3r 

_2 

3 M 
3y 2 39 

3 2 M 

9y 3 90 





TABLE IX.- Continued 


FORTRAN 

Name 


Definition 


D2AY2T 



9y 2 90 


D2AY3P 


DADY1 


DADY2 


DADY3 


DADR 


DADT 


DADP 


DBDR 


DBDT 


DBDP 


9B 

dp 








FORTRAN 

Name 


Definition 


DBDSI 


D2MDSR 


D2MDST 


D2MDSP 


D2UDSR 


D2UDST 


D2UDSP 


D2UY1R 


D2UY2T 


D2UY3P 



YPEMDU 






TABLE IX.- Concluded 


FORTRAN 

Name 

Definition 

YP2MDU 

- i u ~ — 


2 M 9Y 2 

YP3MDU 

Y —LA — C ” — 
3 M 9Y 

O 


DMUDRP 

3l 


9r ' 


2 

DUDSRP 

9 A* 


9 ipdr ' 


2 

DUY1RP 

_ 9 H 


9Y 1 9r> 


2 

DUY2RP 

0-M 


9Y 2 9r ' 


2 

DUY3RP 

9 £ 


9Y 9r ' 

RD0TP 

r ' 

DEL2S 

v 2 ^ 


2 

YD (11) 

V £ 

P = 

F 2 


H cosa 


Subroutine INPUT 


Description.- The purpose of the INPUT routine is to bring 
into core the user's data and then convert it into a useable for- 
mat. (See flow chart in Figure 11.) Most of the initialization 


78 





is performed here. Since all the information referring to the in- 
put section is already presented in Section II, it is sufficient 
to present only the flow chart and the listing of the subroutine. 

Dictionary of major FORTRAN names. - Table X contains a dic- 
tionary of major FORTRAN names. 


TABLE X 


FORTRAN 

Name 

Data 

Type 

Definition 

V 

Numeric array 

Coefficients used to calculate 
N x (r) 

W 

Numeric array 

Coefficients used to calculate 
N x (r) 

D 

Numeric array 

Dummy input array 

XNAME1 

Name list 

Description in Section II 

XNAME3 

Name list 

Description in Section II 

XNAME5 

Name list 

Description in Section II 

H 


Scale size of electron density 
distribution 

HPRIME 


Point of maximum electron den- 
sity gradient 

ZX 


r 

— (r - r ) 
r 

ENF 


N p Electron density in F 

region 

ENXR 


V r > 

ENX 


N Electron density in exo- 

sphere 

EN 


Total electron density, 

N + N = N 
X F 


l 

i 

/ 
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TABLE X.- Concluded 


FORTRAN 

Name 

Data Type 

Definitions 

PKDELN 


AN 0 = PKFRAC * N (PKFRAC is 
in percent) 

PLT 


Array used to store plotting 
title information i 

YO 

Array 

Contains value of differen- 
tial equations evaluated ] 

NUAR 

i 

1 


Number of differential equa- 
tions evaluated 


Subroutine Output 

Description . - Final calculations are made in OUTPUT prior to 
printing results. 


(1) DOPPLER SHIFT - DSHIFT 


Af = -(2.424067E-4) 


A 


I & 

n s<P dh p 


(2) VALIDITY CRITERION - VALCRIT 

2 _ / „ \2 


v c = K [(i) + 7 (fe) + ] 


2 -i 1/2 


where 


K = 


K 1 fy 2 .10 6 


K. = — = 2 . 0943933E-5 
l c 

c = speed of light in km 
f = signal frequency in megacycles/sec. 
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Figure 11B 


Flow Chart of Subroutine INPUT 




/ COMPUTE \ 


V 




J 


/nuar is the'X 

( NUMBER OF ) 
V. EQUATIONS J 


TERM< 0 
N I Y 


COMPUTE 
INTERMEDIATE 
TERMS NEEDED 
BY FIELD 


/fieldS 


COMPUTE 
EMUS = U‘ 


3 


COMPUTE 

INITIAL VALUES 
OF DIFFERENTIAL 
EQUATIONS, STORE 
IN ARRAY YO 

INUAR = 10 


THE VALUE OF 
I TERM IS NEGATIVE! 


3 


( RETURN ) 

GO ON 
,TO NEXT CASE 


) 



EMUS < 0 



N 

Y 

H 

EMUZ 

= (EMUS) 1//2 



THE VALUE OF 
EMUS IS NEGATIVE 


3 


C RETURN ) 

(go ON TO \ 
yNEXT CASE J 


Figure 11C. - Flow Chart of Subroutine INPUT 
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In the flow chart convention: 


DMyDR = 

dr 


Y0(1) = r 


DMUDT = fSJ- 
ou 


RCT = rcos0 


DMUDP = 


.a* 


EMUS = n 


(3) TOTAL POWER LOSS IN DB - PL - 


P = 10 log [fj, PR0PT EXP 


where 




o n cosa 


YD < U > = P F = 2 


fj, cosa 


Y0(11) = P F dh p 


PR0PT = x - 

T 47rr M x 10 1 
o 


(4) GROUP DELAY IN MILLISECONDS - GRODEL 


GR0DEL 


h 

= I I / 1 . f x 1 . 0E6 & \ 

c J 0 ( M df ) 


where f is the frequency in MHz and dh = y cos a ds. 

P 
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In the flow chart convention: 


YD(7) = G t 
L 


+ f • 1.0E6 


af 


Y0(7) 



= G 


L 


GR0DEL = 


G 


L 


c 


where 


G = group path length in km 
c = speed of light in km/sec 
= 3 x 10 2 

Also, see the flow chart. Figure 12. 

Dictionary of major FORTRAN names .- Table XI contains a 
dictionary of major FORTRAN names, subroutine OUTPUT. 


TABLE XI 


FORTRAN 

Name 

Data 

Type 

Definition 

THETA 


Colatitude of wave front 
in degrees 

PHI 


Longitude of wave front 
in degrees 

Y0SQR 


[ Y0 ( 4 ) ] 2 + [Y0 (5) ] 2 + [ Y0 ( 6 ) ] 2 

DSHIFT 


Doppler shift, Af. 

C2 


4 . 24737E-6/F, used to calcu- 
late dispersion 


86 






FORTRAN 

Name 


TABLE XI.- Concluded 


FORTRAN 

Name 

Data 

Type 

Definition 

VALCRIT 


Validity criterion for y 2 

PL 


Total power loss in db, 
exclusive of absorption 

GR0DEL 


Group delay in milliseconds 

XX 

Array 

Array used to store geocen- 
tric radius for plot of 
total ray path 

YY 

Array 

Array used to store colati- 
tude for plot of total ray 
path, theta in degrees 

LL 


Number of points in the XX 
and YY arrays 500 points 
maximum 

JFLAG 


Program indicator, used to 
indicate reason OUTPUT was 
called 

JUMP 


Program indicator, used to 
indicate that an end of case 
condition exists 

N 


Program indicator, used as a 
line counter 


Subroutine COLL 


Description. - The collision frequency model has the followinq 
functional form: 


where v' is computed as shown below. (Also, see flow chart. Fig- 
ure 13.) The collision frequency profile as a function of 
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SUBROUTINE 
OUTPUT j 


CONVERT 
6 and4> 

TO DEGREES 


COMPUTE 

DOPPLER SHIFT (D SHIFT) 
AND 

VALIDITY CRITERION (VALCRIT) 
FOR L 2 


JUMP = 1 


PLOT INTERVAL > . 

HAS BEEN REACHED; | J FLAG = 10 

Y ! N 


N = 0 
Y I N 


STORE PLOTTING DATA 

LL = LL + '1 
XX (LL) = YO ( 1 ) 

YY (LL) = THETA 


LL 2. 500 

N I Y 


JUMP = 

1 I 2 


Figure 12A. - Flow Chart of Subroutine OUTPUT 



















Figure 12C. - Flow Chart of Subroutine OUTPUT 
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^RETURN^ — [ 


N PLOT = 0 
Y I N 


/ CALL \ 

/ CALL \ 
< CALCO 4 > 

t 

RESET COUNTERS 
LL = 0 
LB = 0 
JUMP = 1 

I 


RETURN 


TO GENERATE RECTANGULAR PLOT CALL 
TO GENERATE POLAR PLOT CALL 


CALCO 5 
CALCO 4 


BOTH PLOTS ARE DESCRIBED IN THE FOLLOWING SECTION 


Figure 12D 


Flow Chart of Subroutine OUTPUT 



altitude consists of three parts that are smoothly joined with 
the aid of a curve fitting program. 

For 6378 < r < 6478 km: 

v' = 12.03527 - 0.07392 x 

where 

x = (r-6378) km. 

For 6478 r 6853 km: 

6 

v ' = ^ [a. + C (d, <p) b.J f.(x) 

i=l 


where 

f^x) = 1 

f 2 ( x ) = x 

f 3 ( x ) = x 2 
f 4 (x) = X 3 

f (x) = cos(0.0157 x) 

O 

f„(x) = sin(0.0157 x) 
b 

a x = 5.0562 
a = -3.7482 x 10~ 

dt 

a = 1.3864 x 10~ 4 

O 

a . = -1.4777 x 10‘ 

4 

a_ = -0.48192 

5 

a„ = -0.27021 

6 


; x = (r -6478) km 


; b = 0.032512 
; b 2 = -0.8847 x 10 
; b = 0.8541 x 10 _ ‘ 

O 

; b 4 = -1.5422 x 10 

; b r = 0.01470 

5 

; b = 0.65037 
b 


92 


where 


2 2 

c (0> 0) = Cj0 + c 2 d + c 3 + (d^Q + d 2 0 + dg) cos <j) 


c = -0.35818 
c 2 = 1.1250 
c_ = -0.88344 

U 


; d 1 = -0.17828 
; d 2 = 0.55997 
; d 3 = 0.56028 


0 is the colatitude and <j) is the longitude in degrees. 
<() = 0 corresponds to local noon. For r ^ 6853 km, 

v • = 2.3653 - 0.0030266 x • (0.3195 - 0.0000536 x) c(d, <f>) 


where 


x = (r - 6853) 

Also, see flow chart of subroutine COLL, Figure 13. 

Subroutine FIELD 

Description . - An idealized dipole model is used for the 
magnetic field of the Earth. The magnetic field equations which 
define gyrofrequency , fjj, and the angle between the magnetic 
field and the wave normal, , are 

f H = °Xl (r ) [1 + 3 °° 39 ] 1/2 


where a = 6378 km, the radius of the Earth and r and 0 are the 
geocentric radius and colatitude of the ray position 


c 


11 


t - — B • 10~ 6 = 0.9 

27rm o 


where B 0 = 0.3142 Gauss is the magnetic field on the surface 
of the Earth at the equator and e and m are the charge and mass 
of an electron. 
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\ 


x = : 



Ficrure 13. - Flow Char 



of Subroutine COLL 















cos ^ = —r 


cos0 + YgSinfl 


^ Y l • [1 + 3cos0] 


1/2 


i=l 


2Y 2 cos 0 - Yj^sine 2 ,1/2 

8to * = l2Y 2 co S e-Y 1 sta9| l 1 - °° S *J 


The derivatives of f H and are: 




d(f, 


H 

d0 


= -2.7 


sin© cos0 


£l + 3cos 2 0j 1 ^ 2 


d(f H ) 

d</> 


= 0 


Let 


T ! ■ E Y L 

L = 1 

T g = Y 1 cos0 + | Y 2 sin0 

Tg = cos 2 e + J sin 2 0 
T 4 = - Y 1 sin0 + Y 2 cos0 
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Then 


Also , 


d( co sift) 
dr 


= 0 


d(cosift) 

d0 


T + 3 

" T sin0cos0 

1 + 

4 4 

T 3 


(T 1 > 1/2 • (T 3 ) 1/2 


d( co sift) 
d c p 


= 0 


d(cosift) 

dY 

1 



r y • t i 

1 9 

COS0 - 

T 1 J 


(T ^ 172 (T 3 ) 1/2 


d(cos;ft) 

dY„ 


2 Bin9 - 

r Y 2* T 2 i 

t i 


L 1 J 


(T 1 )1/2 (T 3 )1/2 


d(cosift) 

dY„ 


-Y co sift 

O 


d (cosift) _ 
dY 2 d0 


Y 

d(cos^ft) 2 


4cos0 


d9 T 1 (Tl ) 1/2 ,T 3 ) 1/2 [l + 3co S 2 9 ] 2 


d (sinift) _ 1 

dY 2 d0 simp 


rnof . d (cosift) d(cosift) d (cosift) 

COS! ^dY 2 d0 d0 dY 2 


2 

sin ift 


see flow chart, Figure 14. 
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Dictionary of major FORTRAN names .- Table XII contains a 
dictionary of major FORTRAN names, subroutine FIELD. 


TABLE XII 


FORTRAN 


Name 

Definition 

FH 

f R , gyrofrequency in MHz 

DFHDR 

d ( f H ) 


dr 

DFHDQ 

d(f H> 


d0 

DFHDP 

d V 


d<p 

COSPSI 

cos(ip) 

SP2 

, 2, 

1 - cos (^>) 

DCPDR 

d ( cos #) o 

dr 

DCPDT 

d(cos^) 


d0 

DCPDP 

d(cos^) Q 

d<£ 

DCPDY1 

d(cosip) 

®i 






TABLE XII.- Concluded 


FORTRAN 


Name 

Definition 

DCPDY2 

d(cos^) 


dY 2 

DCPDY3 

d(cos^) 


dY 


3 

D2CY2T 

d(cos^) 

d(Y 2 )d( 0 ) 

D2SY2T 

d(sin^) 

dY o d0 


2 


Subroutine DENSE 


Description . - The electron density is given by 
N(r ,6,<f>) = N^(r, 0) • N 2 (r, 0, <p) 

where r is the geocentric radius, 6 is the colatitude, and 
<j) is the magnetic longitude. 

Ni(r,9) represents the background density in the magnetic 
meridional plane and N 2 (r, 6 ,<|>) represents the ionization of the 
field-aligned irregularities in both the meridional and azimuthal 
planes . 


The background ionization is assumed to have the form: 


N, (r , 0) = m + N 
F x 

where 

N p = electron density in the F-region 
N x = electron density in the exosphere 

N F - N ma X * Exp t 1 /2{l-W- e - W }] 

N max is assume d to be 2.7E5 
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and 


h-h 


W = 


max 


H t 


max 


= altitude of the peak electron density N 


max 


H_ = the scale height 
F 


max 


350 km 


H = 50 km 
F 

The electron density in the exosphere is given by: 


N x = N r x [ 1 . 659E3 • 0 + 1 . 662E4] 


where 


where 


6 


is 


the colatitude in radians. 



1/2 


Z 


x 


is the geopotential altitude given by: 


Z = 
x 


r 

o 

r 


< r - r o> 


where r 0 is the reference radial distance and is equal to 
6878 in the model. 


The symbol i refers to each ion present in the exosphere. 

Hi and 3i refer to the scale height and the fractional density 
of the i-th ion at the reference level r 0 . 

A three-ion gas model consisting of oxygen (O ) , helium (He ) , 
and hydrogen (H"*") ions is assumed. The various parameters in the 
density equation have the following values: 


ion^ 

h 

H i 

0+ 

.9788 

66.546 

H+ 

.0016 

1056.3 

He + 

.0196 

265.98 


The technique suggested by Swayze is used to join the exo- 
spheric profile smoothly with that of the F-layer between 350 and 
1000 km. The expression N x is modulated by a factor 


100 


where h is the altitude in km. 


The model for the field aligned ionization irregularities 
is given by: 


AN 


N 2 (r,0)= 1 + 


K) 


• EXP 


Az 


■•(t) 


N 0 is the peak ionization enhancement and H 0 is the scale 
size of the irregularity at the base of the field line 

t_ _ meridional width at (r, 0) 

t meridional width at (r ,0 ) 
o ' o o' 

sin 3 0 (4-3 sin 2 0 q ) 1 

sin 3 0 o (4-3sin 2 0) 1/2 


Since H p = 50 km 


N = N • exp 
F max 


dN 


[;(-»-*) 


^-■v K) K) 


dN p 

dr 

d N F 

de 


= N„ X 10 2 (e W -l) 


= 0 


N x = N r [K l 9 + K 2 ] 


where K 1 and K 2 are constants and 0 is the colatitude. 


where 


But, 


If N x 


for r 


r 3 


N = 
r 


"1 1/2 


2 >* 


exp 


L i=l 


x 

H. 


Z = ~ (r-r ) 
x r o 


dN 
. x 

dr ' 


(K 2 0 + K 2 ) 


2N 


' 3 

E 

i=l 


- exp 
1 


X 

H. 

l 


dZ 

> 

dr 


_o 

2 


<V + K 2 > 


dZ 

X 

dr 

dN 
x_ 

dr 2N 

i 

dN 

= K.N . 
d0 1 r 

is modulated by the factor 


_o 

2 


E r exp 


i=l 


H. 

l 



2 I 


/ r-7378 \ 

exp 

V 500 t 


H. 

l 


< 7378 km. 


dN 

x 

dr 

dF r 

dr 


d(N -F ) 
v x r 

dr 


dN 


dF 


x r 

— - • F + N — 

dr r x dr 


F . 

r \ 125000/ 
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Similarly, 



d(N 


•F ) 
r' 


60 


= F • 
r 


dN 
x 

6Q 


= K..N • F 
1 r r 


DIS refers to the normal distance of the ray-position (r,0,<t>) 
from a field-line of colatitude X. 


where 



SOA refers to the geometrical path length in Earth radii of any arbitrary ray- 
position (r, 0) from its initial position (r 0 , d Q ) along the same magnetic field line 
of the colatitude \ . 


SOA = 


JE. 

2 

sin X 


/4 2 x 1 / 2 / 4 2 \ 1/ ^ 2 

- sin Gj cosd - I ~ - sin G q I cos 6^ 


+ |log 


(*-"A 


\l/2 


+ COS0 


(I - sl " 2 ®0 ) 


1/2 


+ COS0 


Also, see flow chart. Figure 15. 

Dictionary of major FORTRAN names .- Table XIII contains a 
dictionary of major FORTRAN names, subroutine DENSE. 
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Figure 15A. - Flow Chart of Subroutine DENSE 
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TABLE XIII 


FORTRAN 

Name 

Data 

Type 

Definition 

ENF 


N F - Electron density in the 
F-region 

ENX 


N x - Electron density in the 
exosphere 

ZX 


r 

— (r - r 0 ) - geopotential 
r altitude 

ENXR 


Z V EXP (- 57 j 

DZXDR 


d < Z x> 

dr 

DENFDR 


opy 



dr 

DENXDR 


d < N x> 

dr 

DENXDT 


d < N x» 

d0 

R 

r 

Geocenter radius of wave front, 
in km 

TH 

0 

Colatitude of wave front, in 
radians 

AMBDA 

X 

Colatitude of point where the 
field line intersects the 
Earth's surface, in radians 
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TABLE XIII.- Concluded 


FORTRAN 

Name 

Data 

Type 

Definition 

ZM5 

t 

meridional width 

modula- 

% 

at (r,6) 

meridional width < 

at (r 0 ,e 0 ) 

| tion factor 
v if variable 
| width duct 
is considered 

H 


Scale size of the field line 
at (r , 0 , <j>) 

DIS 


Distance perpendicular from 
(r,9,<t>) to field line 

DD 


Plotting array used to store 
distance from field line, DIS 

S0A 


Plotting array 
distance along 

used to store 
field line 


Subroutine CALC05 

Description . - This subroutine does the plotting for tne ray 
tracing program. All plotting operations are performed by this 
routine . 

Two types of plot are generated by this program. The first 
is a rectangular plot of distance from the field line vs. dis- 
tance along the field line from near end. The polar plot shows 
the total ray path with respect to the Earth. All distances are 
in kilometers. 

Multiple entry points allow the user to obtain either one or 
both of the plots. The user also has the option of overlaying 
successive plots. 

This subroutine has entry points CALC04 and CALC05 in sub- 
routine OUTPUT. CALC01, CALC02 and CALC03 are called by INPUT. 

This subroutine in turn calls subroutine AXIS, SYMBOL, NUM- 
BER, LINE and PLOT. These are standard calcomp plotting sub- 
routines and are described in Bulletin #170-C, Programming for 
Calcomp Digital Incremental Plotters, by California Computer Pro- 
ducts, Inc. 

The subroutine PRAM calculates the adjusted minimum and delta 
required by the line and axis subroutines. 
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Restrictions of CALC04 

Titling - Because only one set of axes is drawn, when using 
the overlay option, titling information for the overlayed plots 
is not printed. The user should know beforehand the content of 
the plot. 

Also, see flow chart. Figure 16. 

Dictionary of major FORTRAN Names .- Table XIV contains a table 
of major FORTRAN names, subroutine CALC05. 


TABLE XIV 


FORTRAN 

Name 


Data 

Type 


Definition 


PLDAT 

X 

Y 

LL 

DATE 

PLT 

XMAXO 

XMINO 

YMAXO 

YMINO 

XMAXl 

XMIN1 

YMAXl 

YMAXl 


Numeric array 

Numeric array 
Numeric array 
Scalar 


Alphanumeric 

array 

Numeric array 


Scalar 


Scalar 


Scalar 


Scalar 


Working storage for plot 
subroutine 

X coordinates to be plotted 

Y coordinates to be plotted 

Number of points in X and Y 
arrays 

20 characters reserved for 
data 

Numeric values of labling 
information 

Maximum and minimum values 
of the X-axis for the rec- 
tangular plot of the dis- 
tance from the field line 
vs. the distance along the 
field line 

Maximum and minimum values 
of the Y-axis for the rec- 
tangular plot of the dis- 
tance from the field line 
vs. the distance along the 
field line 

Maximum and minimum values 
of the X-axis for the polar 
plot of ray path 

Maximum and minimum values 
of the Y-axis for the polar 
plot of ray path 
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Figure 16A. - Flow Chart of Subroutine CALCO 
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Subroutine FORCE 


Description.- Subroutine FORCE solves the equation. 


2 2 2 
Y 1 + Y 2 t Y 3 



In theory this relationship holds all along the ray path but in 
fact it does not. We make it hold by scaling Y 2 and Y 3 down 

proportionately. This process is called "Normalization" of the 
Y vector. FORCE is called by MAIN only. 


Also see flow chart. Figure 17. 



Figure 17. - Flow Chart of Subroutine FORCE 
Subroutine POLAR 

Description . - The expression R = Ejr/E— , a complex number, 
is the polarization term. Ey and Ey are_the components of the 
electric vector_of the wave along axes (x,y) in the wave front, 
y being in and x perpendicular to the plane containing the mag- 
netic field of the Earth. 
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Also, see flow chart. Figure 18. 

Mathematical Subroutines 

CSTl and TOR .- CST1 evaluates the sine and cosine integrals 
Si (x) and Ci (x) . It is called only in MAIN and in turn calls 
TOR. TOR evaluates n! and is called only by CSTl. 

Minv . - MINV solves three sets of three simultaneous equations 
using a Gaussian metnod with pivot selection. It is called by 
POWERL and calls no other subroutine. 

SMARK and MARK .- Subroutines SMARK and MARK constitute the 
integration package. MARK is a MAP language subroutine originally 
written by the Jet Propulsion Laboratory. The purpose of SMARK 
is to allow an interface with the main program written in FOR- 
TRAN IV. 

Listings of the five mathematical subroutines are appended 
in Appendix B. A complete description of the integration package 
is found in Appendix A. 
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Figure 18 


Flow Chart of Subroutine POLAR 
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APPENDIX A 


DESCRIPTION OF THE INTEGRATION PACKAGE 
Introduction 

The integration package is composed of two subroutines, SMARK 
and MARK, both written in the MAP assembly language. Subroutine 
SMARK serves as a connecting link between the FORTRAN IV monitor 
and subroutine MARK. Subroutine MARK as originally written was 
designed to operate in a non-FORTRAN environment. The integration 
of the differential equations is performed by MARK. 

Subroutine SMARK 

This subroutine is buffer between MARK, a differential equa- 
tion solving routine written in MAP language, and a main program 
written in FORTRAN. SMARK allows a program written in FORTRAN IV 
to use most of the features of MARK. 

Usage of SMARK.- A knowledge of MARK is assumed. 

Calling sequences: 

CALL SMARK, KIND, N.HBANK, NRTN , NTRG 
EUBAR, ELBAR, HMAXT , HMINT, YCLOW 
LV1 , TVl , 

LV2 , TV2 , 

LV3 , TV3 , etc., 

up to 10 triggers 


where 

KIND = type of integration 
0 = fixed AM integration 
2 = RK integration 

4 = Variable Adams-Moulton integration 

N = actual number of differential equations 

HBANK = location of a bank of storage used by MARK 
equivalent to HBANK-3 
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NRTN = return indicator from SMARK 

1 = EOS 

2 = DERI 

3 = DER2 

4 = trigger return 

5 = error return 

NTRG = return indicator from SMARK which designates which 
trigger has been activated. 1 for first trigger, 

2 for second, etc. 

EUBAR, ELBAR, HMAXT , HMINT , YCLOW same as in MARK 

LVl = location of variable being tested 

TV1 = location of desired value of the variable tested 

CALL TRA14 returns control to SMARK and causes a TRA 1,4 return 
to MARK. 

CALL TRA24 returns control to SMARK and causes a TRA 2,4 return 
to MARK. 

CALL ON (N) turns trigger N on. 

CALL OFF (N) turns trigger N off. 

The order of the differences to be carried in the Adams-Moulton in- 
tegration must be stored in HBAND(l), the nominal step size in 
HBAND ( 4 ) , the maximum number of equations allowable in HBANK ( 5 ) , 
and all independent and dependent variables initialized before 
calling SMARK. The double precision part of the independent 
variable HBANK (7) is set to zero by SMARK. 


Subroutine MARK 

MARK is a closed subroutine designed to solve the first n 
of a set, N, of first-order differential equations utilizing 
Adams-Moulton open or open and closed formula types. A Runge- 
Kutta fourth-order integrator is used as a starting routine to 
generate backward differences initially. Provision is made for 
interrupting the integration process at specific values of either 
the independent or the dependent variables. The order of dif- 
ferences (m) used in the Adams-Moulton mode is less than or equal 
to nine (9) (m < 9) . 
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Restrictions . - The following restrictions to subroutine MARK 
apply: 

(1) MARK will not integrate backwards in the independent 
variable. The nominal step size, H, must be positive. 
Changes in H must be accomplished by the use of a 
"doubling" or "halving" procedure in MARK that will 
double (set H = 2H) or halve (set H = 0.5H) the inte- 
gration step size. 

(2) Underflow and overflow are not checked internally. 

(3) The user must provide the necessary interruption sub- 
routines, an auxiliary program to evaluate the n 
first-order derivatives, and a bank of storage for in- 
ternal calculation. 

(4) This is an FAP program and is not FORTRAN compatible. 

Integration technique .- The following integration technique 
applies : 

(1) MARK permits the user to solve the N differential 
equations by one of three options: 

a. Runge-Kutta fourth order. 

b. Adams-Moulton with a fixed step size, H, and the 
ability to alter H by the doubling and/or halving 
procedure using Runge-Kutta to initially generate 
backward differences. This applies to either a 
predictor or a predictor with q corrections (open 
or open and closed type formulas) . 

c. Adams-Moulton as mentioned in b. using an automatic 
variable step size control. Halving and doubling 
are controlled automatically. The correction formu- 
la is applied only once. 

(2) Both the independent and the dependent variables are 
automatically carried internally in partial double 
precision to control round-off error locally. The user, 
however, will recognize the variables only as single 
precision quantities. However, the user may carry the 
independent variable in full double precision by option. 
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Usage of MARK .- The following usage of MARK applies: 
Calling Sequence: 

CALL MARK or TSX $MARK , 4 

PZE HBANK , P, EOS 

PZE DERI, <J>, DER2 

ERROR RETURN 

Pfx Bl, , Y1 

PZE Z1 

Pfx B2 , , Y2 

PZE Z2 


Pfx BJ , , YJ 
PZE ZJ 


PZE O 


where the symbols are defined as follows: 


HBANK - The location of a bank of storage to be described 
below. 


0 - The independent variable is carried in partial 

precision. 

P - double precision (single precision to the user) . 

1 - The independent variable is carried in full double 

precision. 

EOS - The location of a user "end of step" routine. This 
routine must terminate with a TRA 1,4 command. It 
is used to evaluate variables that are needed only 
after a full integration step is completed. 
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DERI - The location of the entry to the user's derivative 
routine that carries out all calculations that in- 
volve the independent variable. This routine must 
terminate with a TRA 1,4 command. 

DER2 - The location of the entry to that portion of the 
user's derivative routine that carries out all 
calculations that do not involve the independent 
variable but are required to evaluate the deriva- 
tives. 

A simple example of the use of DERI, DER2 follows: 

Suppose we are to solve 


= 

dx 


ax 


+ by 


then 

DERI ax 2 

Thus , the DERI entry calculates the extra term involving the in- 
dependent variable x. This provides a saving of real machine 
time, particularly during the Runge-Kutta phase of integration, 
and also saves machine time when the closed type formula is used 
with Adams -Moulton integration. 

0 - Adams-Moulton integration with fixed step size 
<p = 2 - Runge-Kutta integration only 

4 - Adams-Moulton using automatic variable step size * 
control 

The pairs of locations in the calling sequence are specified as: 

Pfx BJ, ,YJ 

PZE ZJ are defined as "triggers" 

These triggers are the linkage control to the user's inter- 
ruption subroutines. The triggers state that control is transfer- 
red to location BJ when the contents of location YJ are equal to 
the contents of location ZJ. Thus BJ is the location of a user's 
interruption subroutine, YJ is the location of a variable being 
checked, and ZJ is the location that contains the desired value 
for YJ. 

Triggers are separated into two classes: Independent vari- 

able triggers called T-stops and dependent variable triggers 
called Y-stops. 


121 



( 1 ) 


Independent variable triggers called T-stops. These 
triggers interrupt on values of the independent variable 
of integration. All T-stops must have Y J = 0 . That is, 
they must have the following format in the calling se- 
quence: 

Pfx BJ 

PZE ZJ 

The logic used to execute T-stops is as follows: 

Let tg^, t S 2 , tgj t s k be a set of values of the 

independent variable for which interruptions are desired. 

MARK sets t m = Min [t sl , t s2 , t s },] • Integration continues 

normally until the independent variable reaches the condition: 


The step size is set = (t^ + y -t m ) and integration is carried to 
t m where all the values of the variables including derivatives 
and end of step values are calculated and control is then trans- 
ferred to the user's interruption routine, all values are reset 
to station t^.g and the next t m is determined. If no other t m 
exists within this step, integration continues. Thus, interrup- 
tion routines for all t m within a given step are executed before 
integration continues. There is no limitation on the number of 
T-stops permitted (except for machine size, of course) . 

(2) Dependent variable triggers called Y-stops. These trig- 
gers are interrogated at the beginning of an integration step and 
a value. 


i . = y 


l 


is calculated and saved for each of the j Y-stops. 
of the integration step the difference. 


At the end 


r . = y 


n + 1 " y j 


is calculated and the algebraic sign of r j is compared to 1^ 
If 


sgn 1^ ? sgn r_. 
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then the condition y = yj has occurred within the integration step 
and a linear interpolation search procedure is executed to deter- 
mine the value of the independent variable, t, such that y = y-i . 
When the At calculated by the search procedure is such that J 

At < 6 

u 


-26 

2 Max H, t for P = 0 

rj + 1 

-42 

2 Max H, t for P = 1 

tj + 1 

then convergence to tj is assured. At this point all values of 
the dependent variable including their respective derivatives and 
any end-of-step calculations are determined and control for the 
corresponding Y-stop is returned to the user's interruption rou- 
tine. If more than one Y-stop trigger occurs within an integra- 
tion step, then the triggers are executed in the order of the 
smallest value of the independent variable determined for the 
respective Y-stops. Thus, the order of execution is determined 
by the independent variable. After all Y-stops within an inte- 
gration step have been determined and executed, the conditions at 
station t^ + i are restored for all dependent variables and their 
derivatives and end-of-step calculations, if any. Integration 
then continues normally. 

Up to and including 50 dependent variable triggers are per- 
mitted. However, this number may be altered by changing the sym- 
bolic card "OMAR EQU 50" in the symbolic program deck to the de- . 
sired number. 

It remains to define Pfx of the trigger pair. This is uti- 
lized to permit the user to render triggers "active" or "inactive". 
Active means that a trigger is to be interrogated and executed if 
necessary. Inactive means that the trigger is to be ignored. 

Thus, if 


where 



PZE = trigger is active 

Pfx = 

MZE = trigger is inactive. 

The interruption routines provided by the user must terminate 
with either a TRA 1,4 command or a TRA 2,4 command. 


TRA 1,4 


is used when the interruption does not constitute 
a discontinuity in any of the calculations. 
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TRA 2,4 


is used when a discontinuity exists. Under this 
condition a "restart" procedure is initiated by 
MARK by continuing beyond the discontinuity point 
using Runge-Kutte until a sufficient number of 
backward differences are determined to switch to 
Adams-Moulton integration. 

Comments on triggers .- The following are comments on triggers: 

(1) There is no limitation on how many times a trigger may 
be executed. 

(2) Care must be exercised in updating the ZJ of triggers. 

If the ZJ are not updated after a trigger returns con- 
trol to the user, a machine loop will result, since 
MARK will continue to return control to the user, a 
machine loop will result, since MARK will continue to 
return control to the user's respective interrupt rou- 
tine on the basis of the current Z J. Thus, a trigger 
must either be updated or rendered inactive to looping. 

(3) In all cases where more than one trigger is to be execu- 
ted at a single point (tj) the triggers will be executed 
in order of their ascending appearance in calling sequence 

(4) Control is returned to the error return of the calling 

sequence whenever t m < (t^ - S u ) or when the number 

of Y-stcps exceeds 50. 

The entire list of triggers must be terminated with PZE 0. 

This is the end of the calling sequence for MARK. 

The bank of storage specified by the location HBANK is as 
follows : 

PZE m 
PZE NH 
PZE ND 

HBANK 

DEC H 
PZE N, , n 
DEC 

DEC t 2 
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DEC 

DEC y 2 


DEC y n 


DE c y N 

DEC y^ 
DEC y 2 


DEC 


I 



DEC y^j 


BSS 


3N + 2N(m + 1) for 4 . 0,2 


where 


5N + 3N (m + 2) for <J> 4 


m - order of differences to be carried in the Adams -Moulton 
mode. 
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m ><: 9 (fixed point in the address portion of the word) 
for <J> = 0. 

m ^ 8 for 4> = 4 . 

number of times to sequentially halve the step size in 
the Adams-Moulton mode (fixed point in the address 
portion of the word) . 


NOTE 

NH takes precedence over ND and 
doubling is not executed until the 
number of times to halve is com- 
pleted. If these numbers are in- 
troduced initially in the HBANK, 
the procedure is commenced auto- 
matically when conversion from 
Runge-Kutta to Adams-Moulton is 
completed. NH and ND are ignored 
when using the automatic variable 
step size mode. NH and ND may be 
set by dependent variable or inde- 
pendent variable interruption rou- 
tines in the Adams-Moulton fixed 
mode. Control is returned to the 
user anytime through an interruption 
routine. The number of times halving 
and/or doubling have/has been com- 
pleted is available in the decrement 
portion of NH and/or ND. If addition- 
al halving and/or doubling requests 
are entered in the address portions 
of NH and/or ND before a preceding 
request is completed, the sum of the 
additional request and those remaining 
uncompleted will be executed. 

nominal step size (floating point) . 

total number of 1st order differential equations 
(fixed point) . 

total number of the first n 1st order differential 
equations to be integrated by MARK. n ^ N (fixed point) . 


NOTE 

H and N must not be altered unless a 
restart procedure is executed after 


the inital entry to MARK. Item n may 
be altered after the inital entry to 
MARK through an interruption routine. 

If n is increased, MARK restarts. Care 
should be exercised in setting the ini- 
tial conditions corresponding to the 
additional equations to be integrated. 
If n is decreased, MARK continues nor- 
mally integrating the new n set of dif- 
ferential equations. 


t 1 = single precision value of the independent variable 
in floating point. 

double precision value of the independent variable 
t 2 = in floating point. This must be zero initially if 
P = 0 (single precision) . 


*1 

to 


r 




V. 


values of the N differential equations 
for the dependent variables. The initial 
or starting values must be predetermined 
and set by the user (floating point) . 



to 


I 



values of the derivatives of the dependent 
variables calculated and stored by the user's 
derivative routine (DERI, DER2) . An initial 
pass is executed through DERI, DER2 , and EOS 
by MARK before the integration process is 
commenced (floating point) . 


MARK entry points .- Provision is made through entry points 
to MARK to transmit certain information to MARK or to render 
certain information available to the user that is stored inter- 
nally in MARK: 


HC By using the command 

CLA* $HC 

the user has direct access to the current step size 
being used in the integration process. This is not 
necessarily the nominal step size, H, introduced by 
the user in the HBANK (floating point) . 

NI By using the command 

STO* $NI 
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the user informs MARK that he desires i corrections 
to be performed on the predictor formula used in 
the Adams-Moulton fixed mode of integration. See 
Mathematical Description in this Appendix for des- 
cription of the predictor-corrector formulas being 
used. In the automatic step size control mode i 
is automatically 1, and MARK ignores NI. Thus 1 
correction is made for each prediction in this 
mode (fixed point) . 

TGLO By using the command 

CLA* $TGLO 

The user has direct access to the most recent 
t n + i calculated, where t n + ]_ represents the 
value of the independent variable at the end 
of an integration step (floating point) . 

Y The command 

CLA* $Y 

gives the user access to the location of the de- 
pendent variables (single precision) in the HBANK. 
This appears as L(Y), where index register 1 
set to n and counted down renders all the variables 
to the user (floating point) . 

YDOT The command 

CLA* $YDOT 

y(2) performs the same function as Y for the deriva- 
tives of the dependent variables (floating point) . 

The command 

CLA* $ Y ( 2 ) 

renders the location of the double precision part 
of the dependent variables available to the user 
YO (floating point) . 

The commands 


CLA* $ YO 
CLA* $ YO ( 2 ) 
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render the locations of the single and double pre- 
cision values of the dependent variables at t^ 
available to the user. Item t^ represents the 
value of the independent variable at the beginning 
of an integration step (floating point) . 

The following symbols refer to entry points used for the auto- 
matic step size mode. 

EUBAR The command 

STO* $ EUBAR 

stores E for use in automatic error control (float- 
ing point) . 

ELBAR The command 

STO* $ ELBAR 

stores E in floating point for use with AEC . 

HMAXT 


STO* $ HMAXT 


stores maximum allowable H for AEC (automatic 
error control) (floating point) . 


HMINT 


STO* $ HMINT 


stores minimum allowable H for AEC in floating 
point. 


YCLOW 


STO* $YCLOW 

stores Y for AEC in floating point. 
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RGERR 


CLA* $ RGERR 


permits access to the maximum + q for the user 
in floating point. 


NOTE 

EUBAR through YCLOW 
are consecutive loca- 
tions in MARK. 

Space required .- MARK requires 3453g = 1835qo storage loca- 
tions. No COMMON is required. The user must supply 

5N + 7 + 2N(m + 1) 

storage locations for 4> = 0,2 or 7N + 7 + 3N (m + 2) for <f> = 4 . 

N = maximum number of differential equations; m = order of dif- 
ferences to be carried in the Adams-Moulton mode. 4> = 0,2 is 
for Runge-Kutta integration or for Adams-Moulton integration in 
the fixed mode. Also, whatever storage is required for the user's 
derivative box and trigger control must be supplied. 

Timing information .- MARK will do approximately 40 integra- 
tion intervals per second. (This time was obtained from solving 
a set of 14 first-order differential equations.) 

Checkout . - MARK has been checked out rather extensively using 
a variety of programs at the Jet Propulsion Laboratory. These 
programs include the JPL tracking program, a low thrust traject- 
ory program, and a program of a general nature that solves a sys- 
tem of differential equations starting with five equations, re- 
peating these five and adding sets of five with repetition until 
a maximum of thirty equations have been reached and integrated. 


Mathematical Description 

The classical Runge-Kutta fourth-order equations .- Let the 
system of equations to be solved be in the form, 


y i ' f i “• V-- 


y n ) 3 = 1 > 2 > 


N 
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of^^at VS r^ ^ t n and fj „ be the derivative 

- at t - tr. Let h be the step size of J {±!e independent vari- 


yj 


able t. Then, 


K = h f.(t , y. ) 

1 J J) y i,n 


k 2 = h + 1/2 h - y j,, + V 2 > 

K 3 = ht A + 1/2h - jj,, + V 2) 

K 4 = h V‘t? + At ’ y i,„ + K 3> 


"j, „ + 1 




. The Ada ms-Moulto n predictor-corrector equations 
be defined as above. Then , 


Let 


Yj 'Yj 


Y T n + 1 = y i ». + h(a o V ° + a i ^ + • • • + a V m ) y! 

J,rj + 1 j,T 7 o 1 m ' J j 

(open type) 

where V is a backward difference operator operating on y' and 

j f n 

V°y = y 

* 3f T) 

The predictor coefficients a are* 

m 


a =1.0 
o 

a 1 = 0.5 

a £ = 0.416666666 

a 3 = 0. 375 

a 4 = 0.348611111 

a g = 0.329861111 

a = 0.315591936 
6 
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a = 0.304224539 
7 


a = 0.294868003 
8 

a = 0. 2870754484 
9 


' p 


f,(t , y ) 

] V } 


j = 1, 


N 


y i,v + i = y i.i 


+ h(b V° 

' O 



+ 


, , P 

•- +b m V y j » T7 + 1 


where A is defined as above, 1 is the first corrector application, 
and the corrector coefficients bm are: 


b = 1.0 
o 

b = - 0.5 

b = - 0.0833333333 

b 3 = - 0. 0416666666 

b = - 0.0263888888 
4 


b p = - 0. 01875 

5 

b = - 0. 0142691795 

6 

b = - 0. 0113673950 
7 

b = - 0. 0093565362 

b = - 0. 0078925542 
9 


NOTE 

t) , = 3, _ 3. 

m + 1 m + 1 m 


continuing 



= y. + h(b V° + b V 1 + . 

+ 1 J j,v ° 1 


, , 1 
b V ) y! 
m JtV 


+ 1 


y i,n 


(i + i) _ 


+ 1 


\n + i 


(i) 


+ h cr e 


(i) 
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where 


m 

<T= Z b ; e <‘> ' (1-1) 

l = o m y j.>j + 1 

i is the i correction on the predictor formula. 

The formula for interpolation to interrupt on a dependent 
variable in the Adams-Moulton mode . - The following formula ap- 
plies: 


q j 



j 


where n = 


t 

JL 



>0 


and 


3 


(M-l)(M-2)... Qz-j) 

0 + 1 )! 


j = 1, . . . ,m 


c. = b. + 
J 3 


q.b. . 
i 3-i 


i=0 


j= 1, . . . ,m 


b. 

3 


corrector coefficients described in 2 above. 



m 


y l,ix~ y t,r\ + 1 + 1 + 


m 

^ i) i — 1> • • • > n 
3=1 


Figure A-l describes the configuration. 
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Figure A-l.- Plot of y versus t 


The formula for interpolation to halve the step size (H) 





t = t - nf h n = 1, 2, . . . ; SL = 1/2, 1/3, . . . 

V 

t - n ft h - t 

_Q H = _n£ 

M h 

Let £ = 1/2, then y = - 1/2 where n represents the absolute value 
of the subscript of t in Figure A-2 . 
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Figure A-2 . - Plot of y' versus t 


In the program: 


q o = m7 


(i + M) and q, = q K M ™ ~ *) 

4 k+l q k (fi + k + 1) ( k + 1) 


where k is the absolute value of the subscript of t in Figure A-2. 

e . Automatic step-size control.- The following denotes automatic 
step-size control: 


E = MAX 
77+I 


; D. = MAX y. 


m + 1 


hV + 1 




The expression E n + ^ represents the maximum error in any of 
the dependent variables in the final iterate Yj f ri + 1 ^ ue to trun- 
cation error in the step from t^ to + p Tne user, through the 
entry points, supplies MARK with a set of values to be described 
as follows: 

(1) E - upper bound on the truncation error E^ + 

(2) E - lower bound on the truncation error E^ + 

(3) h max “ maximum allowable value of the step size. 

(4) h m j_ n “ minimum allowable value of the step size. 

(5) y - a constant used to prevent unnecessary reduction 
whenever lyj,n + ll i- s small. £ > 0. 

The step size, h, is doubled, left alone, or halved depending 
on the following inequalities: 


(1) If E^ + i 4 E for m successive steps, the step size, h, 
is set = 2h. 

(2) If E < Eji + < E, the step size, h, is left alone. 

(3) If Ej| + ^ > I, the step size, h, is set = 1/2 h. 

in 1 

The program preserves all the conditions y j , Vj at t n 

and integrates to + p 

If (1) holds, then MARK sets up the doubling procedure and 
integrates m + 1 more steps checking that (1) holds at each step. 
If (1) holds, the doubling procedure is completed and h = 2h. 


If (2) holds, integration continues normally. 

m "j - i 

If (3) holds, then MARK restores Yj,n r n , 
executes the end of step box to restore those values at ^ 


It 


Finally, the halving procedure is executed and h = 1/2 h. 


Thus, it is never necessary on the basis of error control to re- 
start the integration procedure in the Runge-Kutta mode. Item 
E is approximately equivalent to specifying the number of signifi- 
cant figures to preserve locally throughout the integration. Item 
E should normally range from 10“^ to 10“3. y_ may be determined 
by the user but should probably range from lO - ^ to 1. 
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APPENDIX C 


CHECKLIST OF ELEMENTS IN COMMON AND SUBROUTINE 

IN WHICH USED 



M 

A 

I 

N 

0 

I u 

N T 

P 

U 

T 

F 

I 

P E 

U L 

T D 

P c 

F 0 D P A 

0 W E 0 L 

R E N L C 

C R S A 0 

E L E R 4 

/DATA/ 

X 

X 

X X 

X X X X 

RMAX 

X 

X 



RMIN 

X 

X 



TMAX 

X 

X 



TMIN 

X 

X 



PMAX 

X 

X 



PMIN 

X 

X 



PRNT 

X 

X 



RELERR 

Presently not used 


A0 


X 



B0 


X 



R0 


X 

X 

X 

THETA0 


X 

X 

X 

PHI0 


X 

X 

X 

PLUS 

X 

X 


X 

NP0WER 

X 

X 

X 

X 
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0 

I u 

M N 

A P 

I U 

N T 

T 

p 

u 

T 

F F 

I 

E 

L 

D 

p 

0 

0 w 

R E 

c 

E 

R 

L 

D 

E 

N 

s 

E 

c 

P A 

0 L 

L C 

A 0 

R 4 

NPL0T 


X 

X 






J 

X 

X 







NUAR 

X 

X 







N1 

X 








EN 

X 

X 

X 



X 

X 


DNDR 

X 






X 


DNDT 

X 






X 


DNDP 

X 






X 


EMU 

X 

X 

X 


X 




RTYSQR 

X 



X 

X 




F 

X 

X 

X 





X 

F2 

X 

X 

X 



X 



Cl 

X 

X 




X 



FH 

X 

X 


X 


X 


X 

C0SPSI 

X 

X 


X 


X 


X 

STNPSI 

X 

X 


X 


X 


X 

DFHDR 

X 



X 


X 



DCPDT 

X 



X 


X 



DCPDY1 

X 



X 


X 




290 


o P c 

iuffodPa 
mntioweol 

APPERENLC 

I UULCRSAO 

N T T D E L E R 4 


DCPDY2 

X 


X 

X 



DCPDY3 

X 


X 

X 



SP2 

X 


X 




EMUINT 

X X 

X 





EMUS 

X X 

X 


X 



N 

X 

X 





GNU 

X 

X 





MUFLAG 

X 



X 



NTEST 

Presently not used 





/HBANK1/ 

X X 

X 

X 

X X 

X 

X 

M0RDER 

Presently not used 





N0HALF 

Presently not used 





N0D0UB 

Presently not used 





HBANK 

X X 






N0EQ 

Presently not used 





FINVP 

Presently not used 





FINVPI 

Presently not used 





Y0 

X X 

X 

X 

X X 

X 

X 

YD 

X 

X 


X 
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M 

A 

I 

N 

0 

I U F 

N T I 

P P E 

U U 

T T 

F 

L 

D 

p 

0 D 

OWE 
REN 
C R S 

E L E 

c 

P A 

0 L 

L C 

A 0 

R 4 

MA 

Presently not used 




/ CSCR/ 

X 

X X 

X 

X 

X 

c 

X 

X 

X 

X 


RCT 

X 

X 


X 


S 

X 

X 

X 

X 


RST 

X 



X 


Z 

X 




X 

EM 

X 

X 


X 


TERM 

X 

X 


X 


TERM2 

X 

X 


X 


RM0D 


X 



X 

RARG 


X 



X 

Y06 

Presently not used 




/CONST/ 

X 

X 




0RDER 

X 

X 




EUBAR 

X 

X 




ELBAR 

X 

X 




YCL0W 

X 

X 




HMAXT 

X 

X 
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N 

T 

T 

D 

E L 

HMINT 
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